text:
# Generation of polarization-entangled photon pairs in a cascade of two type-I crystals pumped by femtosecond pulses

## I INTRODUCTION

Photon pairs generated in the process of spontaneous parametric down-conversion (SPDC) have been an effective and convenient source of two-particle entangled states used in tests on the foundations of quantum mechanics as well as application to quantum information technologies such as quantum cryptography and quantum teleportation. In such applications, a pulsed source of entangled photon pairs is particularly useful because the times of emission are known to the users. For example, in the entanglement-based quantum key distribution protocol, it is convenient to share entangled photons between sender and receiver in a train of short optical pulses synchronized with the common timing clock pulses. Such a pulsed source allows the sender and receiver to know the arrival times of the photons within the duration of the pump pulse and to time stamp each received bit readily so that they can sift their raw data to generate a shared key from their public discussion. A pulsed source is also crucial in a certain class of experiments that require the use of several photon pairs at a time, such as quantum teleportation, entanglement distillation, and the generation of multiphoton entangled states. In these experiments, time-synchronized entangled photon pairs must be available. It has also been pointed out that primitive elements of a quantum computer can be constructed from the combined use of several entangled photons and Bell-state measurement.

A great deal of effort has been devoted to developing a pulsed source of entangled photon pairs. Over the past decade, femtosecond-pulse-pumped SPDC has been extensively studied by several groups. It has been shown that the entangled states generated in type-II SPDC pumped by ultrashort pulses shows disappointingly low visibility of the quantum interference, and narrow- bandwidth filters are required to increase the visibility at the expense of the photon flux. Theoretical studies revealed that the down-converted photons are entangled simultaneously in polarization and space-time, or equivalently, wave number-frequency because of the significant effects of dispersion and birefringence in the crystal. As a result of this multiple-entanglement, it is possible to distinguish, with some degree of certainty, which of the interfering pathways occurs by the measurement that distinguishes the space-time component of the state, for example, the measurement of the arrival time of the photons at the detector. This distinguishing space-time information is enough to seriously degrade the visibility of two-photon polarization interference. The presence of this distinguishing space-time information stems from the fact that the down-converted photons are emitted spontaneously, but the times of emission are known to be within the duration of the pump pulse. The resultant two-photon state is well localized in space-time, and provide undesirable timing capability. The situation is serious in the ultrashort-pulse-pumped SPDC, but is absent in the ordinary cw-pumped SPDC. Therefore, it is vital in the ultrashort-pulse-pumped SPDC to ensure that the space-time component of the state contained no distinguishing “which-path” information for the photons in order to observe the entanglement and quantum interference in the polarization degree of freedom. Various interferometric techniques have been proposed to eliminate the entanglement in unnecessary degrees of freedom and to recover the quantum interference.

On the other hand, cw-pumped SPDC exhibits a polarization-entanglement with high-visibility quantum interference capability. In particular, Kwiat et al. realized maximally polarization-entangled photons using noncollinear SPDC with type-II phase matching. More recently, they demonstrated that noncollinear SPDC using two spatially separate type-I nonlinear crystals pumped by a cw laser exhibits high-visibility quantum interference. The latter source is particularly convenient since the desired polarization-entangled states are produced directly out of the nonlinear crystal with unprecedented brightness and stability and without the need for critical optical alignment. Their result attracted great attention for the experimentalists. Recently, Kim et al. showed that polarization-entangled photons from two spatially separated type-I nonlinear crystals pumped by femtosecond laser pulses exhibit high-visibility interference. In their experiment, collinear degenerate type-I SPDC and a beam splitter were used to create a polarization-entangled state. The state prepared after the beam splitter was, however, not considered to be entangled without amplitude postselection. Only when one considers those postselected events (half of the total events) in which photons traveled to different output ports, one can observe quantum interference. Later, this problem was solved by using collinear nondegenerate type-I SPDC and a dichroic beam splitter. It was shown that mostly maximally entangled states showing high-visibility (92%) quantum interference were obtained without amplitude and spectral postselection just before detection. However, the available flux of the entangled photon pairs was still very limited ($``$100 s-1).

In this paper, we apply femtosecond-pulse pumping to the second Kwiat scheme with type-I phase-matching arrangement. We shall demonstrate that highly entangled polarization states with sufficiently large flux are obtained when an appropriate temporal compensation is given in the pump pulses. The results will be shown to be consistent with the theory which incorporates the space-time dependence of interfering two-photon amplitudes, if the effects of dispersion and birefringence in the two-crystals are appropriately taken into account.

## II EXPERIMENT

Our experimental setup is schematically shown in Fig. 1. Two adjacent, thin, type-I crystals (BBO) whose optic axes are horizontally (H) and vertically (V) oriented, respectively, are pumped by 45\- polarized femtosecond pulses at 266 nm. The pump pulses were third harmonic of a mode-locked output of Ti/Sapphire laser, whose approximate average power was 150 mW and repetition rate was 82 MHz. Due to type-I coupling, H-polarized photon pairs at 532 nm are generated by the V-polarization component of the pump field in the first crystals, and V-polarized photon pairs are generated by the H-polarization component of the pump field in the second crystal. These two possible down-conversion processes equally likely occurs and are coherent with one another. The SPDC was performed under a degenerate and quasi-collinear condition. The signal and idler photons making an angle of 3 with respect to the pump laser beam and having the same wavelength around 532 nm were observed through the identical interference filters, centered at 532 nm, placed in front of the detectors. The polarization correlations were measured using polarization analyzers, each consisting of a rotatable half- wave plate (HWP) and a quarter-wave plate (QWP) (for 532 nm) followed by a polarizing beam splitter (PBS). After passing through adjustable irises, the photons were collected using 60-mm-focal-length lenses, and directed onto the detectors. The detectors (PMT) were photomultipliers (HAMAMATSU H7421-40) placed at $``$1.5 m from the crystal, with efficiencies of $``$40% at 532 nm and dark count rates of the order of 80 s-1. The photodetection area was about 5 mm in diameter. The outputs of the detectors were recorded using a time interval analyzer (YOKOGAWA TA-520), and pulse pairs received within a time window of 7 ns were counted as coincident.

To obtain a truly polarization-entangled state, care must be taken to disentangle the polarization degree of freedom from any other degrees of freedom, that is, to factorize the total state into product of the polarization- entangled state and those describing other degrees of freedom. This is equivalent to saying that effective polarization entanglement requires the suppression of any distinguishing information in the other degrees of freedom that can provide potential information about “which-polarization” the emitted pair have. In our case, since which crystal the origin of each pair is and their polarization are intrinsically correlated, distinguishing “which-crystal” information must be also eliminated. Accordingly, to make emitted spatial modes for a given pair indistinguishable for the two crystals, we overlap the down-conversion light cones spatially by using very thin crystals each having $``$130 $`\mu `$m-thickness. In addition, it is important to eliminate distinguishing space-time information inherent in the two-photon states produced in the ultrashort-pulse-pumped SPDC. Figure 2 schematically illustrates what may happen when a 45\- polarized femtosecond pulse incidents on two cascaded type-I BBO crystals, which is estimated roughly from the optical characteristics of the BBO. For simplicity, we consider here the degenerate collinear SPDC. In the first crystal, H-polarized down-converted photon wavepackets are advanced at $``$74 fs relative to the V-polarized pump pulse due to group velocity dispersion, while the H-polarized pump pulse is delayed at $``$61 fs relative to the V-polarized pump pulse due to birefringence in the BBO. In the second crystal, V-polarized down-converted photon wavepackets are advanced at $``$74 fs relative to the H-polarized pump pulse due to group velocity dispersion, while the H-polarized down-converted photon wavepackets are advanced at $``$242 fs relative to the H-polarized pump pulse due to an advance given in the first crystal ($``$135 fs) and dispersion in the second crystal ($``$107 fs). Consequently, after the crystals, space-time components of the two-photon state associated with the polarization states $`|HH`$ and $`|VV`$ are expected to be temporally displaced by $``$168 fs, where $`|H`$ ($`|V`$) means a single photon linearly polarized along a horizontal (vertical) axis and the first (second) letter corresponds to the signal (idler). Note that the distance between the two BBO crystals is not an important factor in this consideration. Since down-converted photon wavepackets are expected to have widths comparable to that of pump field ($``$150 fs), it is suggested that the space-time components associated with $`|HH`$ and $`|VV`$ do not overlap in space-time. As a result, “which-polarization” information may be available, in principle, from the arrival time of the photons at the detector.

To eliminate this distinguishing space-time information, a polarization dependent optical delay line for the 266-nm pump was inserted before the crystals, which is denoted as the pre-compensator in Fig. 1. It consists of quartz plates, whose optic axes are oriented either vertically or horizontally, and a Bereck-type tiltable polarization compensator. The combination of the quartz plates of various thicknesses provides a relative delay of $`\left|T\right|`$350 fs between the H\- and V-polarization components of the 266-nm pump field, while the sign is given by the orientation of its optic axis. This delay line can compensate for relative delay between the state associated with $`|HH`$ created in the first crystal relative to the state associated with $`|VV`$ created in the second crystal to overlap them temporally. The Bereck-type polarization compensator is used to adjust the subwavelength delay. As a result, a two-photon Bell state $`|\mathrm{\Phi }^+=(|HH+|VV)/\sqrt{2}`$ can be directly created.

Two kinds of experiments were carried out to confirm whether the target state $`|\mathrm{\Phi }^+`$ was successfully prepared. The first experiment was a quantum state tomography. The polarization density matrices were estimated from 16 kinds of joint projection measurements performed on an ensemble of identically prepared photon pairs. These joint measurements consist of four kinds of projection measurements onto {$`|H,|V,|D,|L`$} on each member of a photon pair, where $`|D=(|H+|V)/\sqrt{2}`$ and $`|L=(|H+i|V)/\sqrt{2}`$. We used a maximum likelihood calculation to estimate the density matrix. We also calculated the concurrence of the state, which is known to give a good measure for the entanglement of a two-qubit system. The second experiment was a conventional two-photon polarization interference experiment. In this experiment, QWPs in mode 1 and mode 2 were set to 0, i.e., their optic axes were vertically oriented, while HWP in mode 1 was set to $`\pm `$22.5 and HWP in mode 2 was rotated. In our setup, this corresponds to the projection measurements onto {$`|L`$ or $`|R`$}$`|\theta `$, where $`|R=(|Hi|V)/\sqrt{2}`$ and $`|\theta =\mathrm{cos}2\theta |H+i\mathrm{sin}2\theta |V`$. The visibility of an observed interference pattern gives another convenient measure for the entanglement.

At first, we used the interference filters having a full-width half maximum (FWHM) bandwidth of 8 nm in front of the detectors. Figure 3 illustrates the estimated polarization density matrices prepared when (a) no relative delay, (b) relative delay of 135 fs, and (c) relative delay of 231 fs were introduced between the H\- and V-polarization components of a 266-nm pump field. In these experiments, maximum coincidence counts and typical accidental coincidence counts were $``$450 cps and $`<`$10 cps, respectively. As shown by these examples, we could successfully prepare the states approximately described by $`\rho (v)=\frac{1v}{2}(|HHHH|+|VVVV|)+v|\mathrm{\Phi }^+\mathrm{\Phi }^+|`$. These states should exhibit interference in coincidence rates for the two detectors in proportion to $`P(\theta )=\frac{1}{4}(1\pm v\mathrm{sin}4\theta )`$. Figure 4 shows the coincidence rates as a function of the HWP angle in mode 2 and the evaluated visibilities associated with the states in Fig. 3(a)-(c). Note that, no background was subtracted in these results, so that the evaluated visibilities give the lower limits of $`v`$. The experiments were repeated for various values of delay $`T`$ introduced in the H\- and V-polarization components of the pump field. The evaluated concurrence C and visibility v of the polarization interference are plotted in Fig. 5 as a function of relative delay $`T`$. This figure clearly indicates that there is a strong correlation between the concurrence and the visibility. This is quite natural because $`C=v`$ should hold for the states $`\rho (v)`$ given by the above form. In this experiment, the maximum values of the concurrence and visibility were 0.95 and 0.92$`\pm `$0.02, respectively. The broken line in Fig. 5 is a Gaussian curve fitted to the measured data, and the FWHM width is 237 fs.

Next, to investigate the effects of interference filters used in front of the detectors, we replaced them with those having a FWHM bandwidth of 40 nm. Figure 6 shows C and v evaluated by the same procedure, except that the backgrounds have been subtracted for both the concurrence and the visibility data in this case. The maximum values of C and v were slightly reduced to 0.92 and 0.90, respectively, while the maximum coincidence counts were increased by a factor of $``$ 6. (Unfortunately, the accidental coincidence increased more than six times, which is probably due to ambient light.) The broken line in Fig. 6 is a Gaussian curve fitted to the measured data, and the FWHM width is $``$205 fs. The width of the fitted curve was slightly reduced in comparison with that for the 8-nm filters, which was due to the spectral filtering effects, as shown later.

We also measured the autocorrelation of the pump field by measuring the degree of polarization, or equivalently, the magnitude of the Stokes parameter $`\left|𝐬\right|=\sqrt{s_1^2+s_2^2+s_3^2}`$ of the pump just after the pre-compensator (and before the crystals), where $`s_1`$, $`s_2`$, and $`s_3`$ are three independent Stokes parameters. Simple analysis shows that $`\left|𝐬\right|`$ is proportional to the temporal overlap of the H\- and V\- polarization components of the pump field between which relative delay $`T`$ was introduced by the pre-compensator, and directly gives the autocorrelation of the pump field. Figure 7 shows the measured $`\left|𝐬\right|`$ as a function of $`T`$. We found that $`\left|𝐬\right|`$ can be fitted with a Gaussian curve having a FWHM width $``$216 fs. Assuming the transform limited Gaussian-shape pulse, the FWHM temporal width of the pump field is estimated to be $``$153 fs, and that of the pump intensity to be $``$108 fs. We see that the width of the autocorrelation of the pump field is comparable to that of the measured concurrence and visibility data shown in Figs. 5 and 6. This result suggests that the width of the measured concurrence and visibility data are determined chiefly by the temporal width of the pump field.

## III THEORETICAL ANALYSIS

Experimental results are analyzed based on the theory given by Grice et al., Keller et al., and Kim et al.. In the following analysis, we will focus our attention on the Gaussian curve fitted with the measured concurrence and visibility data shown in Figs. 5 and 6, which will be denoted by $`V(T)`$. In particular, we will analyze how the widths of these curves are determined and depend on the experimental parameters. To make the following discussion self-contained, it may be helpful to begin by showing the minimum of the previous theories before discussing our problem. Let us consider the degenerate, quasi-collinear SPDC in a single type-I crystal pumped by a field with the central frequency $`\mathrm{\Omega }_p`$ and the FWHM temporal width $`\tau _p^{FWHM}`$. For simplicity we ignore the transverse components of the photon wave vector. Let $`\omega _p[\omega _s,\omega _i]`$ be the frequency and $`k_p(\omega _p)[k_s(\omega _s),k_i(\omega _i)]`$ the wave number for the pump \[signal, idler\]. The perfect phase-matching condition is given by $`\mathrm{\Omega }_s=\mathrm{\Omega }_i=\mathrm{\Omega }_p/2`$ for the frequency and $`K_o(\mathrm{\Omega }_s)+K_o(\mathrm{\Omega }_i)=K_e(\mathrm{\Omega }_p)`$ for the wave number, where $`K_j(\mathrm{\Omega })`$ denotes the wave number for a photon having a frequency $`\mathrm{\Omega }`$ and polarization along the ordinary (j=o) or extraordinary (j=e) optic axis. To simplify the notation, we introduce detuning for the signal, the idler, and the pump from the perfect phase-matching condition by $`\nu _s=`$ $`\omega _s\mathrm{\Omega }_p/2`$, $`\nu _i=`$ $`\omega _i\mathrm{\Omega }_p/2`$, and $`\nu _p=\omega _p\mathrm{\Omega }_p`$, respectively. To first-order in the interaction, the general expression of the two-photon state after the crystal is

$$|\mathrm{\Psi }=\frac{1}{2\pi }𝑑\nu _s𝑑\nu _i\mathrm{\Psi }(\nu _s,\nu _i)a^{}(\nu _s)b^{}(\nu _i)|0,$$
(1)
where $`a^{}(\nu _s)`$ and $`b^{}(\nu _i)`$ are the photon creation operators for the o-polarized signal and idler modes characterized by detuning frequencies $`\nu _s`$ and $`\nu _i`$, respectively, which are defined after the crystal. The function $`\mathrm{\Psi }(\nu _s,\nu _i)`$ is the two-photon amplitude in the frequency domain, and is given by

$$\mathrm{\Psi }(\nu _s,\nu _i)=0\left|a(\nu _s)b(\nu _i)\right|\mathrm{\Psi }=C\alpha (\nu _s+\nu _i)\mathrm{\Theta }(\nu _s,\nu _i).$$
(2)
In this expression, $`\alpha (\nu _p)`$ is the spectral envelop function of the time-dependent classical pump field at the input surface of the crystals, $`\mathrm{\Theta }(\nu _s,\nu _i)`$ is the longitudinal phase-matching function, and C is a constant. Assuming a Gaussian shape for the temporal envelop function of the pump field, $`\alpha (\nu _p)`$ is given by

$$\alpha (\nu _p)=e^{\left(\frac{\nu _p}{\sigma _p}\right)^2},$$
(3)
where $`\sigma _p=\sigma _p^{FWHM}/2\sqrt{\mathrm{ln}2}`$ $`=4\sqrt{\mathrm{ln}2}/\tau _p^{FWHM}`$ is the FWHM bandwidth of the pump field. If z-direction is taken to be the pump direction, the longitudinal phase-matching function $`\mathrm{\Theta }(\nu _s,\nu _i)`$ is given by

$$\mathrm{\Theta }(\nu _s,\nu _i)=\frac{1}{L}_{L/2}^{L/2}𝑑ze^{i\mathrm{\Delta }(\nu _s,\nu _i)z}=sinc[\frac{\mathrm{\Delta }(\nu _s,\nu _i)}{2}L],$$
(4)
where L is the crystal length, $`\mathrm{\Delta }(\nu _s,\nu _i)=k_s(\nu _s+\mathrm{\Omega }_p/2)+k_i(\nu _i+\mathrm{\Omega }_p/2)k_p(\nu _p+\mathrm{\Omega }_p)`$ is the phase mismatch and approximately given by

$$\mathrm{\Delta }(\nu _s,\nu _i)=D_+(\nu _s+\nu _i)+\frac{1}{4}D^{\prime \prime }(\nu _s\nu _i.)^2.$$
(5)
In this expression, $`D_+`$ and $`D^{\prime \prime }`$ are crystal parameters: $`D_+1/u_o(\mathrm{\Omega }_p/2)1/u_e(\mathrm{\Omega }_p)`$ is the group velocity mismatch where $`u_o(\mathrm{\Omega }_p/2)`$ ($`u_e(\mathrm{\Omega }_p)`$) is the group velocity of the o-polarized down-converted photons (the e-polarized pump photon) inside the crystal and $`D^{\prime \prime }d^2K_o/d\mathrm{\Omega }^2|_{\mathrm{\Omega }=\mathrm{\Omega }_p/2}`$ is the group velocity dispersion for the down-converted photons.

Since we are interested in the temporal characteristics of the down-converted photons, we next consider the two-photon state after the crystal in the time domain. Formally, it can be written as

$$|\mathrm{\Psi }=𝑑t_s𝑑t_i\widehat{\mathrm{\Psi }}(t_s,t_i)a^{}(t_s)b^{}(t_i)|0,$$
(6)
where $`a^{}(t_s)`$ and $`b^{}(t_i)`$ are Fourier transforms of $`a^{}(\nu _s)`$ and $`b^{}(\nu _i)`$. For example,

$$a^{}(t_s)=\frac{1}{\sqrt{2\pi }}𝑑\nu _se^{i\left(\nu _s+\mathrm{\Omega }_p/2\right)(t_sL_s/c)}a^{}(\nu _s),$$
(7)
where $`L_s`$ is the optical path length from the output face of the crystal to the observation point. Physically, $`t_s`$ and $`t_i`$ are interpreted as the time of arrival of the idler and signal at the points that are $`L_s`$ and $`L_i`$ away from the output face of the crystal, and $`\widehat{\mathrm{\Psi }}(t_s,t_i)`$ gives information about the space-time motion of the two-photon wavepackets. The two-photon amplitude in the time domain is simply given by the Fourier transform of the two-photon amplitude in the frequency domain and by

$`\widehat{\mathrm{\Psi }}(t_s,t_i)`$ $`=`$ $`0\left|a(t_s)b(t_i)\right|\mathrm{\Psi }`$ (8)
$`=`$ $`{\displaystyle \frac{e^{i\mathrm{\Omega }_pt_s}}{2\pi }}{\displaystyle 𝑑\nu _s𝑑\nu _ie^{i(\nu _st_s+\nu _it_i)}\mathrm{\Psi }(\nu _s,\nu _i)}.`$ (9)
We can further change the variables according to $`\nu _+=`$ $`\nu _s+\nu _i=\nu _p`$, $`\nu _{}=\nu _s\nu _i`$, $`t_+=(t_s+t_i)/2`$, and $`t_{}=(t_st_i)/2`$, and obtain the equivalent Fourier transform relation

$`\widehat{\psi }(t_+,t_{})`$ $`=`$ $`{\displaystyle \frac{e^{i\mathrm{\Omega }_pt_s}}{2\pi }}`$ (11)
$`\times {\displaystyle }{\displaystyle }d\nu _+d\nu _{}e^{i(\nu _+t_++\nu _{}t_{})}\psi (\nu _+,\nu _{}),`$
where $`\psi (\nu _+,\nu _{})=\mathrm{\Psi }(\frac{\nu _++\nu _{}}{2},\frac{\nu _+\nu _{}}{2})`$ and $`\widehat{\psi }(t_+,t_{})=\widehat{\mathrm{\Psi }}(t_++t_{},t_+t_{})`$. Physically, $`t_{}`$ and $`t_+`$ mean half the difference in the arrival time, and the mean arrival time of the idler and signal photons, at the points that are $`L_s`$ and $`L_i`$ away from the output face of the crystal, respectively. The absolute square of the two-photon amplitude $`\left|\widehat{\psi }(t_+,t_{})\right|^2`$ is proportional to the probability of arrival of the idler and signal at the given points for the given mean arrival time $`t_+`$ and half the difference in the arrival time $`t_{}`$. We show schematically the relevant amplitudes in the time domain and the frequency domain in Fig. 8.

Now, let us proceed to provide a phenomenological model that is consistent with our experiment using two crystals. We assume down-conversion processes are equally likely to occur in two crystals whose optic axes are orthogonally oriented and thy are coherent with one another. We represent the overall two-photon state after the crystals as a superposition of the states arising from the SPDC in each crystal with appropriate phase. Then the two-photon state after the crystals in the time domain can be written as

$`|\mathrm{\Psi }`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle 𝑑t_+𝑑t_{}}`$ (14)
$`\times \{\widehat{\psi }_{HH}(t_+,t_{})a_H^{}(t_++t_{})b_H^{}(t_+t_{})`$
$`+e^{i\varphi }\widehat{\psi }_{VV}(t_+,t_{})a_V^{}(t_++t_{})b_V^{}(t_+t_{})\}|0`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle 𝑑t_+𝑑t_{}}`$ (17)
$`\times \{\widehat{\psi }_{HH}(t_+,t_{})|H,t_++t_{}_s|H,t_+t_{}_i`$
$`+e^{i\varphi }\widehat{\psi }_{VV}(t_+,t_{})|V,t_++t_{}_s|V,t_+t_{}_i\}`$
where the first (second) term represents the two-photon state generated in the first (second) crystal, and $`\varphi `$ is the phase difference between the states generated in the two crystals. In Eq. 17 we incorporated a new dichotomic variable in the creation operators and the two-photon amplitudes that specifies the polarization (H or V) of the photon, and wrote $`|H,t_++t_{}_s=a_H^{}(t_++t_{})|0`$ etc., which denotes a single-photon state for the signal (or idler) mode characterized by a definite polarization and time. From the previous discussion, two two-photon amplitudes $`\widehat{\psi }_{HH}(t_+,t_{})`$ and $`\widehat{\psi }_{VV}(t_+,t_{})`$ are expected to be temporally displaced by $`\tau `$ in the $`t_+`$-direction due to dispersion and birefrengence in the crystals. On the other hand, the pre-compensator for the pump before the crystals introduces temporal displacement in the $`t_+`$-direction for these amplitudes by $`T_p`$. Therefore, we postulate that these two amplitudes can be written using a single common amplitude as $`\widehat{\psi }_{HH}(t_+,t_{})=\widehat{\psi }_{VV}(t_++T,t_{})\widehat{\psi }(t_++T,t_{})`$, where $`T=\tau T_p`$. In this case, the two-photon state is now entangled both in polarization and in space-time, in other words, it is doubly entangled. The space-time components $`\widehat{\psi }_{HH}(t_+,t_{})`$ and $`\widehat{\psi }_{VV}(t_+,t_{})`$ provide the distinguishing information for the interfering path of the two-photon amplitudes and degrade polarization entanglement. The full density matrix associated with the state given in Eq. 17 can be written as

$$\rho _{ijkl}(t_1,t_2,t_3,t_4)=_ij,t_2|_si,t_1|\mathrm{\Psi }\mathrm{\Psi }|k,t_3_s|l,t_4_i,$$
(18)
where $`i,j,k,l=H`$ or $`V`$. Since our chief concern is the polarization state of the photon pair, we will calculate the effective density matrix of the two-photon state after the two crystals that are defined in the polarization space alone. It can be calculated by partially tracing the density matrix in Eq. 18 over the time-variable, and is given by

$$\rho _{ijkl}=𝑑t_1𝑑t_2\rho _{ijkl}(t_1,t_2,t_1,t_2).$$
(19)
Explicit calculation gives the effective density matrices of the form

$`\rho `$ $`=`$ $`{\displaystyle \frac{1}{2}}(|HHHH|+|VVVV|`$ (21)
$`+v(T)e^{i\varphi }|HHVV|+v^{}(T)e^{i\varphi }|VVHH|),`$
where we used the normalization condition of the two-photon amplitude:

$$𝑑t_+𝑑t_{}\left|\widehat{\psi }_{jj}(t_+,t_{})\right|^2=1$$
(22)
and notation $`\rho _{ijkl}=ji\left|\rho \right|kl`$ ($`i,j,k,l=H`$ or $`V`$). In Eq. 21, $`v(T)`$ is the convolution of the two-photon amplitudes in the time domain

$`v(T)`$ $`=`$ $`{\displaystyle 𝑑t_+𝑑t_{}\widehat{\psi }_{HH}(t_+,t_{})\widehat{\psi }_{VV}^{}(t_+,t_{})}`$ (23)
$`=`$ $`{\displaystyle 𝑑t_+𝑑t_{}\widehat{\psi }(t_++T,t_{})\widehat{\psi }^{}(t_+,t_{})}.`$ (24)
It should be noted that $`v(T)`$ is, in general, a complex number. However, in the present experiment, we always make $`v(T)e^{i\varphi }`$ a real number $`\left|v(T)\right|`$ by using the Bereck-type polarization compensator. Then, the density matrix given in Eq. 21 is formally the same as what we observed in the experiment, and Eq. 24 indicates that their concurrence is determined by the convolution of the two-photon amplitudes in the time domain.

Based on the above theory, we have made several numerical calculations to elucidate the physics behind the measured $`V(T)`$ curves. We calculated two-photon amplitudes in the frequency and time domains using the experimental parameters: $`D_+570`$ fs/mm, $`D^{\prime \prime }855`$ fs2/mm, $`L=0.13`$ mm, and $`\tau _p^{FWHM}=153`$ fs. Figure 9 shows (a) the (unnormalized) two-photon amplitude $`\psi (\nu _+,\nu _{})`$ and (b) $`\widehat{\psi }(t_+,t_{})`$. We can see that the bandwidth of $`\psi (\nu _+,\nu _{})`$ is much wider in the $`\nu _{}`$-direction than in the $`\nu _+`$-direction, while the temporal width of $`\widehat{\psi }(t_+,t_{})`$ is much wider in the $`t_+`$-direction than in the $`t_{}`$-direction. This is due to the Fourier-transform relations between $`\psi (\nu _+,\nu _{})`$ and $`\widehat{\psi }(t_+,t_{})`$. We can also see that the temporal width of $`\widehat{\psi }(t_+,t_{})`$ in the $`t_+`$-direction is $`\tau _p^{FWHM}`$. Using these results, the magnitude of the convolution $`\left|v(T)\right|`$ is calculated and plotted in Fig. 10 together with the fitted curve of autocorrelation of the pump field shown in Fig. 7. We can see that $`\left|v(T)\right|`$ nearly agrees with the autocorrelation curve of the pump field. This indicates that the temporal shape of $`\widehat{\psi }(t_+,t_{})`$ is approximately Gaussian in the $`t_+`$-direction and their widths are nearly equal to the temporal width of the pump field. This is consistent with the fact that, although there is a small difference in the width, the measured $`V(T)`$ curves nearly agree with the autocorrelation curve of the pump field as in Figs. 5-7.

Now, let us examine in more detail the origin of the different widths of the measured $`V(T)`$ curves in Figs. 5 and 6 in which the different bandwidth spectral filters were used. We will show that this was due to spectral filtering effects. So far, the two-photon amplitudes after the crystal and before the spectral filter have been discussed. They are intrinsic ones in the sense that they are determined only by the pump field and the crystal parameters and irrelevant to the spectral filters. In contrast, what we actually observed in the experiment was the two-photon state after the spectral filters. Hence, it is required to relate the state after the filters to the state before the filters to find the spectral filtering effects. This can be done straightforwardly. Suppose that the linear spectral filters used for the signal and idler have complex frequency responses expressed by $`F_s(\nu _s)`$ and $`F_i(\nu _i)`$, respectively. These functions are assumed to be normalized such that

$$𝑑\nu _s\left|F_s(\nu _s)\right|^2=𝑑\nu _i\left|F_i(\nu _i)\right|^2=1$$
(25)
holds. Then the two-photon amplitude after the filters are simply given by $`\mathrm{\Psi }_F(\nu _s,\nu _i)=\mathrm{\Psi }(\nu _s,\nu _i)F_s(\nu _s)F_i(\nu _i)`$. Since the two-photon amplitude in the time domain is related to that in the frequency domain by the Fourier-transform relations, the associated two-photon amplitude in the time domain is written as

$`\widehat{\mathrm{\Psi }}_F(t_s,t_i)`$ $`=`$ $`{\displaystyle \frac{e^{i\mathrm{\Omega }_pt_s}}{2\pi }}{\displaystyle 𝑑\nu _s𝑑\nu _ie^{i(\nu _st_s+\nu _it_i)}\mathrm{\Psi }_F(\nu _s,\nu _i)}`$ (26)
$`=`$ $`{\displaystyle \frac{1}{\left(2\pi \right)^2}}{\displaystyle 𝑑\tau _s𝑑\tau _i\widehat{\mathrm{\Psi }}(\tau _s,\tau _i)\widehat{\mathrm{\Phi }}(t_s\tau _s,t_i\tau _i)},`$ (27)
where $`\widehat{\mathrm{\Phi }}(t_s,t_i)=\widehat{F}_s(t_s)\widehat{F}_i(t_i)`$ is a normalized two-variable function that incorporates the temporal response of the filters, and $`\widehat{F}_s(t_s)`$ and $`\widehat{F}_i(t_i)`$ are the Fourier transform of $`F_s(\nu _s)`$ and $`F_i(\nu _i)`$, for example,

$$\widehat{F}_s(t_s)=\frac{1}{\sqrt{2\pi }}𝑑\nu _se^{i\left(\nu _s+\mathrm{\Omega }_p/2\right)t_s}F_s(\nu _s).$$
(29)
Equation LABEL:eq17 indicates that the two-photon amplitude after the filters is given by the convolution of the two-photon amplitude before the filters $`\widehat{\mathrm{\Psi }}(t_s,t_i)`$ and the apparatus function for the filters $`\widehat{\mathrm{\Phi }}(t_s,t_i)`$. Because of this relation, the input function $`\widehat{\mathrm{\Psi }}(t_s,t_i)`$ is smoothed by the filter function $`\widehat{\mathrm{\Phi }}(t_s,t_i)`$ to yield broadened output function $`\widehat{\mathrm{\Psi }}_F(t_s,t_i)`$. It may be worth pointing out that Eq. LABEL:eq17 is formally analogous to the connection of phase-space probability distribution with the Wigner functions describing system and measuring apparatus (quantum ruler or filter). In this case, the concurrence and visibility are determined by the convolution of the output function $`\widehat{\psi }_F(t_+,t_{})=\widehat{\mathrm{\Psi }}_F((t_s+t_i)/2,(t_st_i)/2)`$ in contrast to Eq. 24, in which the input function $`\widehat{\psi }(t_+,t_{})`$ appears. Accordingly, their widths depend on the temporal response of the filters, and equivalently, depend on the bandwidth of the filters. In general, the smaller the filter bandwidth is, the larger the width of the $`V(T)`$ curve is, and vice versa.

Now, let us discuss quantitatively the difference in the width of the measured $`V(T)`$ curves associated with the 8-nm and 40-nm filters. If we assume the frequency response functions $`F_s(\nu _s)`$ and $`F_i(\nu _i)`$ to be well approximated by the Gaussian function with FWHM widths $`\mathrm{\Delta }\nu `$, the temporal widths of the associated response functions $`\widehat{F}_s(t_s)`$ and $`\widehat{F}_i(t_i)`$ is in proportion to $`\mathrm{\Delta }\nu ^1`$. In our case, $`\mathrm{\Delta }t_{8nm}104`$ fs for the 8-nm filter and $`\mathrm{\Delta }t_{40nm}20.8`$ fs for the 40-nm filter. On the other hand, from the experimental results shown in Figs. 5 and 6, those consistent with the present theory are roughly estimated to be $`\mathrm{\Delta }t_{8nm}^{}86.0`$ fs and $`\mathrm{\Delta }t_{40nm}^{}17.2`$ fs. Here, we used the fact that the width $`\sigma `$ of the convolution of two Gaussian functions each having width $`\sigma _1`$ and $`\sigma _2`$ is given by $`1/\sigma ^2=1/\sigma _1^2+1/\sigma _2^2`$. Therefore, even if there remain small quantitative discrepancies, it seems reasonable to conclude that there is considerable validity in the above simplified theory.

Let us switch our attention to the spectral characteristics. For the two-photon state given in Eq. 1, the joint spectral intensity is given by $`I(\nu _s,\nu _i)\left|\mathrm{\Psi }(\nu _s,\nu _i)\right|^2`$. Then the spectra $`I_s(\nu _s)`$ and $`I_i(\nu _i)`$ of the signal and the idler can be obtained by tracing $`I(\nu _s,\nu _i)`$ over the unobserved variable, for example,

$$I_s(\nu _s)𝑑\nu _iI(\nu _s,\nu _i).$$
(30)
Figure 11 shows the calculated individual spectra for the down-converted photons just after the crystal. This figure indicates that the FWHM bandwidth of the down-converted photons is $`40`$ nm. Thus, in our experiment using the 8-nm filters in front of the detectors, what we have observed is only a fraction of the photons having the wavelength $`532\pm 4`$ nm that can reach the detectors. On the contrary, we can safely say that almost all the photons produced in the SPDC have been observed in our experiment using the 40-nm filters. In other words, the coincidence count rate ($`2700`$ cps) obtained in the experiment using the 40-nm filters had already reached the upper limit in our configuration. However, since no optimization was made for the transversal spatial mode of the down-converted photons in this experiment, we believe that further improvement in the coincidence count rate may be possible by optimizing the experimental configuration.

While our simplified model seems satisfactory in essence, there still remains a problem which needs to be solved. Within our model, the autocorrelation curve of the pump field should have the smallest width, whereas this is not the case. Although we have not arrived at a conclusion yet, there might be things neglected in this consideration that should actually have been taken into consideration. For example, we neglected the effect higher than the second-order optical nonlinearity in the crystal. There actually might occur higher order effects. In particular, it is likely that self-phase modulation (SPM) due to third-order optical nonlinearity occurs when the crystals are pumped by ultrashort optical pulses having large peak intensities. When the SPM is taken into account, the two-photon state after the crystal is chirped, and their coherence length might be shortened, and the resultant width of the measured $`V(T)`$ curves might be reduced. To obtain a more satisfactory explanation, the effect of SPM in the crystals might be incorporated into the theory.

## IV CONCLUSION

In conclusion, we have generated pulsed polarization-entangled photon pairs by femtosecond-pulse-pumped SPDC in a cascade of two type-I crystals. It was found that highly entangled photon pairs were successfully obtained by using giving an appropriate temporal delay between the orthogonal polarization components for the pump, without the need of spectral post-selection using narrow-bandwidth filters. Theoretical analysis showed that entanglement depends on the convolution of the space-time components of the interfering two-photon amplitudes. The experimental results obtained were well explained if the dispersion and birefringence effects in the two-crystals are taken into account. Our analysis also revealed the effect of the spectral filtering on the magnitude of the entanglement when the interfering two-photon amplitudes has a different space-time dependence.

This work was supported by CREST of JST (Japan Science and Technology Corporation).


format: text

text:
# Convergence of threshold estimates for two-dimensional percolation

## I Introduction

Percolation StaufferAharony92  is one of the most fundamental and widely studied systems in statistical physics. Theoretical studies of percolation models and applications of percolation theory to physical systems have spawned thousands of papers over the last few decades. Even so, there are some substantial gaps in our understanding of percolation. For example, we have at present no exact value for the position $`p_c`$ of the percolation threshold for site percolation on that simplest of two-dimensional lattices, the square lattice. And in three dimensions we have almost no exact results whatsoever. Because of this, numerical methods have played an important role in the study of percolation. In this paper we consider a class of methods for estimating $`p_c`$ for site percolation using finite-size scaling, and show how various estimates of $`p_c`$ in this class scale with varying system size in two dimensions.

The methods studied here for measuring $`p_c`$ are widely used and are all based upon consideration of the crossing probability function $`R_L(p)`$ StaufferAharony92 ; ReynoldsStanleyKlein80 ; Kirkpatrick80 . This function gives the probability that a connected path crosses the system from one boundary segment to another, at site occupation probability $`p`$ and system size or length-scale $`L`$. Some examples of these estimates are:

1. The renormalization-group (RG) fixed-point estimate $`p_{\mathrm{RG}}(L)=p^{}(L)`$, where $`p^{}`$ is the solution to the equation ReynoldsStanleyKlein80 ; YoungStinchcombe75

$$R_L(p)=p.$$
(1)
2. The average value of $`p`$ at which crossing first occurs StaufferAharony92 ; ReynoldsStanleyKlein80 :

$`p_{\mathrm{av}}(L)=p`$ $`=`$ $`{\displaystyle _0^1}pR_L^{}(p)dp`$ (2)
$`=`$ $`1{\displaystyle _0^1}R_L(p)dp,`$
where the last equality follows from integrating by parts. The prime indicates differentiation with respect to $`p`$.
3. The estimate $`p_{R_c}(L)`$ corresponding to the point where $`R_L(p)`$ equals its universal infinite-system value $`R_cR_{\mathrm{}}(p_c)`$ (which is determined by the system shape and boundary conditions) Ziff92 :

$$R_L(p)=R_c.$$
(3)
For a square system, where $`R_c=\frac{1}{2}`$, this estimate $`p_{0.5}(L)`$ corresponds to the median of the distribution $`R^{}(p)`$. A related estimate $`p_{h+v}(L)`$ for rectangular systems is the value of $`p`$ at which the horizontal and vertical crossing probabilities sum to unity YonezawaSakamotoHori89 ; LanglandsPPS92 ; LanglandsPS94 :

$$R_L^{(h)}(p)+R_L^{(v)}(p)=1.$$
(4)
This estimate is identical to $`p_{R_c}(L)=p_{0.5}(L)`$ when the boundary is a perfect square.
4. The estimate $`p_{\mathrm{max}}(L)`$ which is the value of $`p`$ at which $`R_L^{}(p)`$ reaches a maximum (or equivalently, where $`R_L(p)`$ is at its inflection point) ReynoldsStanleyKlein80 :

$$R_L^{\prime \prime }(p)=0.$$
(5)
5. The cell-to-cell RG estimate, which is the point where two systems of different size have the same value of $`R`$ ReynoldsStanleyKlein80 . One possible choice for this estimate, $`p_{\mathrm{c}\mathrm{c}}^{(1)}(L)`$, is the value of $`p`$ at which

$$R_L(p)=R_{L1}(p),$$
(6)
while a second, $`p_{\mathrm{c}\mathrm{c}}^{(2)}(L)`$, is the point at which

$$R_L(p)=R_{L/2}(p).$$
(7)

In order to use these estimates to determine the threshold precisely, we need to know the manner in which they converge to $`p_c`$ as $`L\mathrm{}`$. While it is possible to simulate very large systems for which finite-size effects may be quite small, the statistics for such simulations are still relatively poor because of the small number of samples that can typically be generated. In most cases, better results can be derived by doing more simulations on smaller systems, and this requires that the finite-size behavior is characterized accurately.

It is usually assumed (based upon very general scaling arguments) that all finite-system estimates of the percolation threshold converge to the bulk value $`p_c`$ as

$$p_{\mathrm{est}}(L)p_ccL^{1/\nu }$$
(8)
where $`c`$ is a system-dependent constant, and $`\nu `$ is the exponent governing the correlation length $`\xi `$, such that $`\xi |pp_c|^\nu `$. (For the two-dimensional systems we will be looking at in this paper $`\nu =\frac{4}{3}`$.) Well known exceptions to this behavior are a few highly symmetric, self-dual systems, such as bond percolation on a square lattice with a square boundary, and site percolation on a triangular lattice with a rhomboidal boundary; in both these cases $`R_L(p)`$ is perfectly symmetric about $`p=\frac{1}{2}`$ for all $`L`$ and all the estimates above give $`p_c=\frac{1}{2}`$ exactly. For these systems, one could consider the constant $`c`$ above to be zero.

In Ref. Ziff92, , however, it was argued that for non-self-dual systems with a square boundary, such as site percolation on a square lattice (where because of the non-duality the estimates show finite-size effects), the convergence of most of those estimates is faster than given by Eq. (6). This is an observation of some practical significance, since this particular system (site percolation on a square lattice with a square external boundary) is one of the most commonly studied systems in percolation. Similar arguments also apply to other symmetric two-dimensional crossing problems, such as a system with a rhomboidal boundary, which is commonly used when simulating triangular and honeycomb lattices.

The arguments of Ref. Ziff92,  were based upon the hypothesis that

$$R_L(p)f_0(x)+L^1f_1(x)+\mathrm{}$$
(9)
for large $`L`$, where $`x=(pp_c)L^{1/\nu }`$, $`f_0(x)`$ represents the universal part of $`R`$, and $`f_1(x)`$ represents the first-order correction to the scaling limit. The choice of $`L^1`$ as the leading order of the correction was based on numerical measurements of $`R`$ at $`p_c`$, and can be derived from the assumption that the system is effectively slightly rectangular in shape, because of the different types of boundary conditions applied along the two principal axes Ziff96 . For small $`x`$, it was assumed that

$`f_0(x)`$ $`=`$ $`a_0+a_1x+a_3x^3+\mathrm{}`$ (10)
$`f_1(x)`$ $`=`$ $`b_0+b_1x+b_2x^2+\mathrm{}`$ (11)
where $`a_0=\frac{1}{2}`$ by the symmetry and self-duality of the square boundary. The same symmetry also implies that $`f_0\frac{1}{2}`$ is an odd function in $`x`$ and hence that $`a_n=0`$ for all even $`n>0`$, as above. (To see that $`f_0\frac{1}{2}`$ is odd, note that for the perfectly dual system of bond percolation on a square lattice, $`R\frac{1}{2}`$ is an odd function of $`x`$ for all $`L`$, including $`L=\mathrm{}`$, and by universality systems with other underlying lattices must behave the same way in the scaling limit.) In Ref. Ziff92, , no particular assumption was made about the behavior of $`f_1(x)`$, other than its analyticity about $`x=0`$.

Note that the form of $`f_0(x)`$ in Eq. (10) is not entirely universal, because the independent variable $`x`$ should incorporate a metric factor which depends upon the underlying lattice AharonyHovi94 ; HuLinChen95 . For convenience, however, since we are considering only the one system of site percolation on a square lattice in this paper, we do not include that factor here.

Once the above assumptions, Eqs. (10) and (11) are made, the convergence of the various estimates of $`p_c`$ is straightforward to analyze, and one finds that while the RG estimate does indeed converge according to Eq. (8) (a result that has been verified in many numerical studies), the rest of the estimates above should converge according to the faster behavior

$$p_{\mathrm{est}}(L)p_ccL^{11/\nu },$$
(12)
where the constant $`c`$ varies from estimate to estimate. Simulations reported in Ref. Ziff92,  for systems of size up to $`1024\times 1024`$ sites verified this convergence for the estimate $`p_{0.5}`$ to high accuracy. The estimates $`p_{\mathrm{c}\mathrm{c}}`$, $`p_{\mathrm{max}}`$ and $`p_{\mathrm{av}}`$ were studied in Ref. Ziff92,  using only exact enumeration results for systems of sizes up to $`7\times 7`$ which give polynomials for $`R`$ (see the appendix), and while the behavior of these results was found to be roughly consistent with (12), the uncertainly due to higher-order corrections was large.

Following the publication of Ref. Ziff92, , Aharony and Hovi AharonyHovi94 ; HoviAharony96  argued that the irrelevant scaling variables in the renormalization-group treatment of percolation imply a slower leading-order convergence of $`R_L`$ to its infinite-system value, characterized by an exponent $`\omega `$, whose value was deduced from the Monte-Carlo work of Stauffer Stauffer81  to be about $`\omega =0.85`$. (Note that Aharony and Hovi used $`\theta _1`$ to denote the exponent we call $`\omega `$.) A variety of series expansion results from the early 1980s were also analyzed to give values for this exponent ranging from 0.89 to over 1 AdlerMoshePrivman82 ; AdlerMoshePrivman83 .

The argument given by Aharony and Hovi implies that the leading terms in the expansion of $`R_L(p)`$ should in fact be

$$R_L(p)f_0(x)+L^\omega f_\omega (x)+L^1f_1(x)+\mathrm{}$$
(13)
where

$$f_\omega (x)=c_1x+c_2x^2+\mathrm{}$$
(14)
Hovi and Aharony argued that the constant term $`c_0`$ is zero for a square system, because at $`p_c`$ there are no correction terms for the square-lattice, bond-percolation system, and corrections due to irrelevant variables should be universal. They also argued that $`f_1(x)`$ should be even, so that $`b_1=0`$ in (11), also by symmetry and self-duality. They discussed various consequences of these assumptions, and presented numerical evidence that the term containing the exponent $`\omega `$ is indeed the leading correction term, by showing that the behavior of $`R_L(p)`$ for large $`x`$ (that is, for $`pp_c`$) was better fit with such a term than without it. However, the procedure they used did not allow them to determine the value of $`\omega `$ accurately, because of the numerical difficulty of finding $`R_L(p)`$ precisely for all $`p`$. Furthermore, they did not study the convergence of the other estimates given above.

Recently NewmanZiff00 ; NewmanZiff01 , the present authors have shown that quantities like $`R_L(p)`$ can be studied efficiently for all $`p`$ by first finding the crossing probability $`R_{L,n}`$ in a “microcanonical” system of exactly $`n`$ occupied sites, and then convolving with a binomial distribution to derive results for the corresponding “canonical” system thus:

$$R_L(p)=\underset{n=0}{\overset{N}{}}\left(\genfrac{}{}{0pt}{}{N}{n}\right)p^n(1p)^{Nn}R_{L,n},$$
(15)
where $`N=L^2`$ for site percolation on a square lattice. The microcanonical crossing probability is found using an efficient cluster-joining algorithm employing data structures based on trees and a fast method is employed for checking for percolation on the fly during the progress of the calculation. (While many of the ideas incorporated in this method were put forward previously Hu92 ; GouldTobochnik96 ; deFreitasLucenaRoux99 ; ShchurVasilyev00 ; Tarjan75 ; JanLookmanStauffer83 , this was the first time that all of these components were combined in this way, for the purpose of finding $`R`$ efficiently.) In Ref. NewmanZiff00,  we studied the function corresponding to $`R`$ for the probability of a cluster wrapping around the boundary of a periodic system on a torus. In the present paper, we describe how that method can be implemented for crossing of an open system, and we report results from some large-scale simulations. The results allow us to determine accurately the behavior of all of the estimates above, and to test the theoretical predictions that follow from Eq. (13). As we will see, the appearance of the “irrelevant” term in the scaling of some estimates is confirmed, and a new, more precise value of $`\omega `$ is found.

The outline of the paper is as follows. In Section II we derive the expected scaling behavior of the various estimates of $`p_c`$, assuming the form (13). In Section III we describe our numerical method, and in Section IV we present the results of our simulations. In Section V we give our conclusions.

## II Convergence of estimates

If we assume Eq. (13) to be a correct description of the behavior of $`R_L(p)`$, it is straightforward to deduce the resulting convergence of the various estimates for $`p_c`$. In the following paragraphs we derive the leading correction term for each estimate, or leading two terms when their powers are close to each other.

1. The RG fixed-point: The relevant terms of Eq. (1) are

$$\frac{1}{2}+a_1x+a_3x^3+\mathrm{}=p_c+xL^{1/\nu },$$
(16)
which implies

$$p_{\mathrm{RG}}(L)=p_c+\left[\frac{p_c\frac{1}{2}}{a_1}\frac{a_3(p_c\frac{1}{2})^3}{a_1^4}+\mathrm{}\right]L^{1/\nu }+\mathrm{O}(L^{2/\nu }).$$
(17)
The term in brackets is the value of $`x`$ that is the solution to $`f_0(x)=p_c`$.
2. Average-probability estimate: Eq. (2) gives

$$p_{\mathrm{av}}(L)=1L^{1/\nu }_{x_0}^{x_1}\left[f_0(x)+L^\omega f_\omega (x)+L^1f_1(x)\right]dx,$$
(18)
where $`x_0=p_cL^{1/\nu }`$ and $`x_1=(1p_c)L^{1/\nu }`$, which are the values of $`x`$ at $`p=0`$ and 1 respectively. Noting that, since $`f_0(x)`$ is odd about $`f_0(0)=\frac{1}{2}`$ and approaches 1 as $`x\mathrm{}`$, we have

$$L^{1/\nu }_{x_0}^{x_1}f_0(x)dxL^{1/\nu }x_1=1p_c,$$
(19)
and we then get

$$p_{\mathrm{av}}(L)=p_c+L^{\omega 1/\nu }_{\mathrm{}}^{\mathrm{}}f_\omega (x)dx+L^{11/\nu }_{\mathrm{}}^{\mathrm{}}f_1(x)dx,$$
(20)
where we have extended the limits of the integrals to $`\pm \mathrm{}`$. This result is also implied by Eq. (40) of Ref. HoviAharony96, , for $`n=1`$. The order of the next correction depends upon the higher-order corrections to (13).
3. Median-$`p`$ estimate: Eq. (3) gives

$$\frac{1}{2}+a_1x+b_0L^1+c_1xL^\omega \mathrm{}=\frac{1}{2},$$
(21)
which implies

$$p_{0.5}(L)=p_c\frac{b_0}{a_1}L^{11/\nu }+\mathrm{O}(L^{1\omega 1/\nu }).$$
(22)
4. Maximum estimate: Eq. (5) gives

$$6a_3x+2b_2L^1+(2c_2+6c_3x)L^\omega +\mathrm{}=0,$$
(23)
which implies

$$p_{\mathrm{max}}(L)=p_c\frac{c_2}{3a_3}L^{\omega 1/\nu }\frac{b_2}{3a_3}L^{11/\nu }+\mathrm{O}(L^{2\omega 1/\nu }).$$
(24)
5. Cell-to-cell estimate: Eq. (6) gives

$`a_0+a_1(pp_c)L^{1/\nu }+b_0L^1+c_1(pp_c)L^{1/\nu \omega }+\mathrm{}`$
$`=a_0+a_1(pp_c)(L1)^{1/\nu }+b_0(L1)^1+c_1(pp_c)(L1)^{1/\nu \omega }+\mathrm{}`$ (25)
which implies

$$p_{\mathrm{c}\mathrm{c}}^{(1)}(L)=p_c+\frac{b_0}{a_1}\nu L^{11/\nu }+\mathrm{O}(L^{1\omega 1/\nu }).$$
(26)
Likewise, for $`p_{\mathrm{c}\mathrm{c}}^{(2)}`$ we have

$$p_{\mathrm{c}\mathrm{c}}^{(2)}(L)=p_c+\frac{b_0}{a_1(12^{1/\nu })}L^{11/\nu }+\mathrm{O}(L^{1\omega 1/\nu }).$$
(27)

Thus, the scaling of $`p_{\mathrm{RG}}`$, $`p_{0.5}`$, and $`p_{\mathrm{c}\mathrm{c}}`$ is unaffected to leading order by the presence of the exponent $`\omega `$. However, $`p_{\mathrm{av}}`$ and $`p_{\mathrm{max}}`$ are affected by $`\omega `$ to leading order, scaling as $`L^{\omega 1/\nu }`$, slightly slower than predicted in Ref. Ziff92,  (order $`L^{11/\nu }`$). Note that $`b_1`$ does not enter in any of these results, so the leading scaling does not change if $`b_1`$ is equal to zero, as argued to be the case in Ref. HoviAharony96, .

## III Procedure

We have performed simulations to test the scaling hypotheses above using the algorithm described in Refs. NewmanZiff00,  and NewmanZiff01, . Briefly, sites are occupied one by one in random order starting with an empty lattice. Occupied sites form contiguous clusters, each of which is identified uniquely by the site label of a chosen single site within the cluster, which we call the “root site.” Other (non-root) sites within a cluster possess pointers that point either directly to the root site, or to other sites within the cluster such that by following a succession of pointers one can get from any site to the root. A newly added site is considered to be a cluster of size one, which is its own root site, and bonds are then added between it and any adjacent occupied sites. The clusters to which sites at either end of such a bond belong are identified by following pointers from them to their corresponding root sites, and if the root sites found are different we conclude that two different clusters have been joined by the addition of the bond. We represent this by adding a pointer from the root site of one of the clusters to the root site of the other. Smaller clusters are always made subclusters of larger ones, and all pointers followed are subsequently changed to point directly to the root of their own cluster. The net result is an algorithm that can calculate $`R_{L,n}`$ (and many other observable quantities) for all values of $`n`$ in average running time which is of the order the area of the lattice, or $`\mathrm{O}(L^2)`$ for a square lattice.

In our previous calculations using this algorithm we measured the probability of the existence of a cluster that wraps around the periodic boundary conditions of a toroidal lattice. In this paper we are interested instead in the existence (or not) of a cluster that spans an open system along one given direction. There are (at least) three efficient methods for detecting spanning of this kind, two of which are described in detail in Ref. NewmanZiff01,  and all of which we have used in the present work. In the first method, we use the same pointer-based trick that we used in Ref. NewmanZiff00,  to detect wrapping with periodic boundary conditions, but start out with an $`(L+1)\times (L+1)`$ lattice in which one horizontal row of sites is fixed to be permanently empty and one vertical one is fixed occupied. Occurrence of a wrapping cluster in such a system is then exactly equivalent to occurrence of a spanning cluster in the horizontal direction in an open system with dimensions $`L\times L`$.

In the second method, two complete rows of sites at the top and bottom of an open $`(L+2)\times (L+2)`$ lattice are fixed permanently empty, and two columns of $`L`$ sites on the left and the right sides of the lattice are fixed occupied. The two columns of occupied sites form two initial clusters on the lattice. By following pointers from one site in each of these clusters it is then simple to determine whether the two clusters have been amalgamated by other sites added between them: if they have the same root site they have been amalgamated, otherwise they have not. Performing this check after the addition of each site to the lattice, we can detect when a spanning cluster on the $`L\times L`$ open lattice first appears. The two methods have comparable running times and give compatible results. The second is somewhat simpler to implement.

In the third method, which was used for the majority of the simulations here, we consider an $`L\times L`$ open lattice and keep track of the minimum and maximum $`x`$\- and $`y`$-coordinates for sites in each cluster, updating their values as necessary when clusters are joined. When $`x_{\mathrm{max}}x_{\mathrm{min}}=L1`$ for a cluster, we know that the cluster spans the lattice in the horizontal direction, and similarly for vertical crossing. This method allows one to check for both crossing events simultaneously, and it is also efficient and easy to program. A similar method was used in Ref. Machtaetal96,  for simulations of the Ising model.

Since in the present calculation we are only interested in the existence or not of a system-spanning cluster, we can stop the simulation once a spanning cluster is detected, as spanning must also occur for all higher values of $`n`$. This produces about a 40% saving in running time. Each simulation then produces just a single number, the value of $`n`$ at which a spanning cluster first appears (or two numbers if we check for spanning in both the horizontal and vertical directions). Making a histogram of these values over many runs of the algorithm, we derive an estimate of the probability $`P_{L,n}`$ that the system first percolates when the number of occupied sites reaches $`n`$. This probability is related to the desired function $`R_{L,n}`$ according to $`P_{L,n}=R_{L,n}R_{L,n1}`$, and hence

$$R_{L,n}=\underset{n^{}=0}{\overset{n}{}}P_{L,n^{}}.$$
(28)
Once the $`R_{L,n}`$ is determined, $`R_L(p)`$, the corresponding function in the “canonical” percolation ensemble, is calculated from Eq. (15), with the binomial coefficients for large $`N`$ being calculated by iterative multiplication NewmanZiff01 . The estimates $`p_{\mathrm{RG}}`$, $`p_{0.5}`$, and $`p_{cc}`$ for the percolation threshold are then evaluated directly according to Eqs. (2), (3), and (6). The remaining estimate, $`p_{\mathrm{av}}`$, could be found directly by performing a numerical integral over $`R_L(p)`$, but a better method is to use the following exact formula:

$`p_{\mathrm{av}}`$ $`=`$ $`1{\displaystyle _0^1}R_L(p)dp`$ (29)
$`=`$ $`1{\displaystyle \underset{n=0}{\overset{N}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)R_{L,n}{\displaystyle _0^1}p^n(1p)^{Nn}dp`$
$`=`$ $`1{\displaystyle \frac{1}{N+1}}{\displaystyle \underset{n=0}{\overset{N}{}}}R_{L,n}.`$
Using Eq. (28) this can also be written directly in terms of $`P_{L,n}`$ as

$$p_{\mathrm{av}}=1\frac{1}{N+1}\underset{n=0}{\overset{N}{}}\underset{n^{}=0}{\overset{n}{}}P_{L,n^{}}=\frac{1}{N+1}\underset{n=0}{\overset{N}{}}nP_{L,n},$$
(30)
which means that the canonical average position of the percolation threshold is $`N/(N+1)`$ times the microcanonical average $`(1/N)nP_{L,n}`$ and no convolution is necessary to find its value. Higher moments of the distribution $`R^{}`$ can be found in a similar fashion. For the second moment, for example, we have

$`p^2`$ $`=`$ $`{\displaystyle _0^1}p^2R_L^{}(p)dp`$ (31)
$`=`$ $`1{\displaystyle \frac{2}{(N+1)(N+2)}}{\displaystyle \underset{n=0}{\overset{N}{}}}(n+1)R_{L,n}`$
$`=`$ $`{\displaystyle \frac{1}{(N+1)(N+2)}}{\displaystyle \underset{n=0}{\overset{N}{}}}n(n+1)P_{L,n}.`$
To find where $`R_L^{\prime \prime }(p)=0`$ for the estimate $`p_{\mathrm{max}}`$, we make use of the following result:

$`R_L^{\prime \prime }(p)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{N}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)p^n(1p)^{Nn}R_{L,n}\times `$
$`\left[{\displaystyle \frac{n(n1)}{p^2}}{\displaystyle \frac{2n(Nn)}{p(1p)}}+{\displaystyle \frac{(Nn)(Nn+1)}{(1p)^2}}\right].`$
The above three results, along with Eq. (15), demonstrate further the advantage of calculating $`R_L(p)`$ through the microcanonical $`R_{L,n}`$: quantities like $`R_L(p)`$ and $`R_L^{\prime \prime }(p)`$ can be calculated exactly at all $`p`$, while $`p_{\mathrm{av}}`$ can be found without introducing any error through numerical integration.

In the more familiar binary search method for finding $`R_L(p)`$ StaufferAharony92 , $`p`$ is increased or decreased to narrow the bounds on one’s estimate of position of the percolation point for a given realization of the disorder. This search process is stopped after some number $`m`$ of iterations, typically about 15, giving to a resolution of $`2^m`$ on the estimate for $`p_c`$. Thus, $`R_L(p)`$ is evaluated at only a finite set of points, and this adds some uncertainty to the calculation, beyond the basic statistical error. Given that our microcanonical method is also much faster than binary search due to its efficient cluster-merging and percolation-checking, there seems no reason to use other methods.

## IV Results of simulations

Simulations were carried out for $`L=7`$, 8, 16, 32, 64, 128, 256. The results are given in Table 1, along with exact results from exhaustive enumeration of states for small systems with $`L7`$. The polynomials from which the exact results are derived are listed in the appendix. We also conducted some simulations at $`L=7`$ to compare numerical and exact results, and the agreement was found to be perfect within the statistical accuracy of the simulations. The pseudo-random number generator used for the simulations was the four-tap feedback generator known as R9689 or gfsr4 Ziff98 .

Error analysis for the simulation data indicates that the estimates of $`p_c`$ are accurate to about six figures, and the values of $`R(p_c)`$ are accurate to four or five figures, as indicated the table. We also simulated $`2.4\times 10^7`$ samples for a system of size $`L=512`$, but the statistical accuracy of the results was insufficient to add anything to the present analysis.

Consider the results for the estimate $`p_{\mathrm{max}}`$, whose convergence is non-monotonic. Its value starts below $`p_c`$ for small systems, then goes above $`p_c`$ as lattice size passes though $`L100`$, and presumably converges to $`p_c`$ from above as $`L\mathrm{}`$. Indeed, according to Eq. (24), $`p_{\mathrm{max}}`$ has two correction terms with closely-spaced scaling exponents $`\omega 1/\nu 1.65`$ (using the value of $`\omega `$ from below) and $`11/\nu =1.75`$; it appears that these terms contribute more or less equally in the range of system sizes that we are considering. The observed behavior is consistent with Eq. (24) if $`b_2<0`$, $`c_2>0`$, and the two are are roughly comparable in magnitude. (Note that $`a_3<0`$ because $`f_0^{}(x)`$ is at a maximum at $`x=0`$). It is not possible to fit the exponents of Eq. (24) reliably to these data.

The rest of the estimates are all monotonic, and lead to reasonably straight lines when viewed on a logarithmic plot of $`|p_{\mathrm{est}}p_c|`$ vs. $`L`$, reflecting the leading power-law behavior. To provide a more sensitive representation of our data, we calculate successive slopes between pairs of points for systems of size $`\frac{1}{2}L`$ and $`L`$, giving values $`\left[\mathrm{log}|p_{\mathrm{est}}(L)p_c|\mathrm{log}|p_{\mathrm{est}}(\frac{1}{2}L)p_c|\right]/\mathrm{log}2`$ for the leading exponent for the various estimates $`p_{\mathrm{est}}(L)`$. For these calculations we used the value $`p_c=0.5927462`$ given in Ref. NewmanZiff00, , which is consistent with the data presented here, but of somewhat higher precision than these data would yield.

In Fig. 1 we show the plots of the successive slopes for estimates $`p_{0.5}(L)`$, $`p_{\mathrm{c}\mathrm{c}}^{(2)}(L)`$, and $`p_{\mathrm{av}}(L)`$, as a function of $`L^\omega `$ with $`\omega =0.9`$. According to Eqs. (22) and (27), both of these estimates should converge with a leading exponent of $`11/\nu `$ and a next-order term of order $`11/\nu \omega `$, which implies that the successive slopes should fall on a straight line when plotted as a function of $`L^\omega `$, with intercept of 1.75. This behavior is indeed seen in Fig. 1, with measured intercepts of 1.754 and 1.763 respectively. (Note that agreement is not highly sensitive to the value of $`\omega `$; if the data were plotted as a function of $`L^1`$, the fit to linearity would not be much worse.)

The successive slopes for $`p_{\mathrm{av}}(L)`$ do not fall on as good a straight line as the other estimates, presumably because the exponents $`\omega 1/\nu `$ and $`11/\nu `$ of the two leading terms in the convergence are closely spaced (see Eq. (20)). Extrapolation to $`L=\mathrm{}`$ is still possible however, and we find an intercept at $`1/\nu +\omega =1.65\pm 0.02`$, clearly different from the value of $`1.75`$ for the other estimates, implying $`\omega =0.90\pm 0.02`$, the figure we have used above. The error bars are smaller than the size of the symbols for all $`L`$ except $`L=256`$, as shown in the plot. Note that although $`\omega `$ is used in the abscissa of Fig. 1, its value has little effect on the determination of $`\omega `$ from the intercept of $`p_{\mathrm{av}}(L)`$.

We can also compare the coefficients for the leading behavior of estimates to their predicted values. For example, the predicted value of the leading coefficient for $`p_{0.5}`$, Eq. (22), is $`b_0/a_1=0.423`$, using values of $`b_0=0.322`$ (see below) and $`a_1=0.765`$ Ziff92 ; HoviAharony96 . This compares favorably with the value measured here of $`0.436`$. And for $`p_{\mathrm{c}\mathrm{c}}^{(2)}`$, the coefficient from (27), $`b_0/[a_1(12^{1/\nu })]`$, should have a value of $`1.04`$, which compares favorably with the measured value $`1.02`$. Note that if we take the linear combination of these two estimates (whose finite-size corrections are opposite in sign), $`(p_{0.5}+\alpha p_{\mathrm{c}\mathrm{c}}^{(2)})/(1+\alpha )`$ where $`\alpha =12^{1/\nu }`$, the leading correction terms are predicted to cancel one another and the combination should have leading scaling of $`L^{1\omega 1/\nu }=L^{2.65}`$. And indeed this combination is seen to converge very quickly in our numerical results, with values 0.592698, 0.592739, 0.592745, and 0.592746 for $`L=32`$, 64, 128, and 256 respectively. The plot of these figures vs. $`L^{2.65}`$ given in Fig. 2 shows linear behavior as expected, with an intercept at $`p=0.5927464(5)`$, consistent with the best current figure of $`p_c=0.5927462(1)`$ NewmanZiff00 . Thus, by taking a combination of estimates, we can improve the convergence rate for the open system to the point where it becomes competitive with that of the periodic system, in which the estimate with the best convergence has exponent $`\frac{11}{4}`$ NewmanZiff00 . This cancellation of leading order corrections between the two terms is expected to be universal.

The results for $`p_{\mathrm{RG}}`$ and the standard deviation $`\sigma =\sqrt{(\mathrm{\Delta }p)^2}=\sqrt{p^2p^2}`$ converge to the infinite system value with exponent $`1/\nu =0.75`$ as shown in Fig. 3 where we plot the successive slopes as a function of $`1/L`$. (The convergence of the standard deviation of the distribution is predicted in Ref. HoviAharony96,  and also follows from the equations of Section II). The intercepts are at 0.749 and 0.763 respectively. We also show in Fig. 3 a plot of the type introduced by Stauffer, in which $`p_{\mathrm{RG}}`$ is shown as a function of $`\sigma `$, allowing extrapolation to infinite system size to be carried out without knowledge of the value of $`\nu `$ StaufferAharony92 . This plot shows nearly linear behavior, with the last three points ($`L=64`$, 128, and 256) well fit by the line $`p_{\mathrm{RG}}=0.5927465+0.231512\sigma `$ with $`R^2=0.9999979`$. The intercept is in excellent agreement with the known value of $`p_c`$, although this agreement is perhaps fortuitous, considering the slow convergence of the RG estimate.

The last column of Table 1 gives the crossing probability at $`p_c`$. The considerations in Section II imply that $`R_L(p_c)1/2+b_0/L`$ with no contributions from the irrelevant scaling variable Ziff92 ; Ziff94 ; HoviAharony96 , and indeed an analysis of this data shows good agreement with the behavior $`R_L(p_c)=0.320/L0.44/L^2+\mathrm{}`$, yielding $`b_0=0.320\pm 0.001`$. This is nearly identical to the value $`0.319`$ given in Ref. Ziff92,  (where larger systems, but with lower statistics, were generated) and the value $`0.31\pm 0.01`$ of Ref. HoviAharony96, .

## V Conclusions

We have studied the finite-size scaling of estimates of the percolation threshold $`p_c`$ derived from the crossing probability $`R_L(p)`$ for site percolation on the square lattice. Our numerical results confirm that different estimates converge to $`p_c`$ with a variety of scaling exponents as predicted by the scaling theory developed in Refs. Ziff92, , AharonyHovi94, , and HoviAharony96, . In particular, we have shown that the average threshold estimate $`p_{\mathrm{av}}`$ converges with a nontrivial exponent $`L^{\omega 1/\nu }`$ whose origins lie in the irrelevant variables in the renormalization group treatment of the problem, and our results for this estimate provide us with a direct measurement of that exponent. We find $`\omega =0.90\pm 02`$, somewhat higher than the value of about 0.85 found previously Stauffer81  but consistent with the (wide) bounds set by series studies AdlerMoshePrivman82 ; AdlerMoshePrivman83 . Our result is in good agreement with a renormalization-group result of $`\omega =0.915`$ found by Burkhardt Burkhardt80 ; note1 . More extensive simulations could be done to give $`\omega `$ to higher precision.

Of the estimates considered here, $`p_{\mathrm{av}}`$ is the only one that shows the effect of the irrelevant exponent clearly. The estimates $`p_{0.5}`$ and $`p_{\mathrm{c}\mathrm{c}}`$ are confirmed to converge as $`L^{11/\nu }`$, as proposed previously Ziff92 . The maximum estimate $`p_{\mathrm{max}}`$ is found to exhibit non-monotonic behavior, which can be explained by competition between correction terms with closely spaced exponents $`11/\nu `$ and $`\omega 1/\nu `$.

The numerical results reported here were found using a microcanonical simulation method which allows one to calculate $`R_L(p)`$ easily for any $`p`$ NewmanZiff00 ; NewmanZiff01 . The various estimates can then be found quickly to any desired degree of precision by applying appropriate formulas, Eqs. (17). This method proves particularly advantageous for the estimate $`p_{\mathrm{av}}`$, since this estimate depends on knowing $`R_L(p)`$ for all values of $`p`$, the determination of which by most other methods requires a great deal of work. From the microcanonical data, $`p_{\mathrm{av}}`$ can be found without any calculational bias using Eq. (30).

Having characterized the convergence rates of our various threshold estimates, one can go back to older literature and find many instances where an anomalous rate of convergence was seen but not fully understood or recognized. For example, Reynolds et al. ReynoldsStanleyKlein80  derived an estimate of $`p_c`$ which is equivalent to our $`p_{\mathrm{av}}`$ from numerical data for the one-way crossing probability which they denoted $`R_1`$. Assuming this estimate to scale as Eq. (8), they plotted their results against $`L^{1/\nu }`$ (their Fig. 13). The resulting plot is seen to fall on a nearly vertical line, consistent with higher-order behavior. Fitting their data with the supposed $`L^{1/\nu }`$ scaling, they deduced a best estimate of $`p_c=0.5931`$ in the large system-size limit. If however one assumes instead the $`L^{\omega 1/\nu }`$ scaling derived here, the intercept of their data becomes $`p_c0.5927`$, which is much closer to the current best estimate of this quantity.

Yonezawa et al. YonezawaSakamotoHori89  plotted a quantity essentially equivalent to our $`p_{0.5}`$ as a function of $`L^{1/\nu }`$, and found apparent agreement with this assumed $`L`$-dependence (their Figs. 7 and 8). The expected $`L^{11/\nu }`$ behavior is evidently too weak to be distinguished from $`L^{1/\nu }`$ within the errors on their numerical data. Similarly, Hu et al. HuChenWu96  believed the cell-to-cell estimate $`p_{\mathrm{c}\mathrm{c}}^{(1)}`$ to be insensitive to $`L`$; again, the precision of their work did not allow them to observe the higher-order scaling predicted by Eq. (27).

It should be emphasized that the convergence behavior discussed here is specific to a two-dimensional system with a square or rhomboidal open boundary, with crossing defined as a path from one specified side to its opposite. Because of the symmetry of this system, some terms cancel out, allowing various higher-order corrections to become dominant. For different boundaries (such as rectangular ones), and for higher dimensional systems, the behavior will generally be different. In those cases, most of our estimates will scale as the conventional $`L^{1/\nu }`$, except perhaps the estimates $`p_{\mathrm{c}\mathrm{c}}`$ and $`p_{R_c}`$. (To employ the latter, $`R_c`$ must be known, but we have exact values only for rectangular and conformally related two-dimensional systems Cardy92 .) Study of these other systems is a subject for future research.

Another approach to measuring $`p_c`$ is to use periodic rather than open boundary conditions. A partially periodic system in two dimensions is a cylinder, and crossing in this system was studied in Ref. HoviAharony96, . The fully periodic rectangle is a torus, and the criterion of crossing is replaced by criteria involving the different topologically distinct ways in which clusters can wrap around the boundaries Pinson94 . (Some authors HuLinChen95 ; Hu96 ; Machtaetal96  have also considered the percolation criterion in which a cluster has the full dimension of the lattice along at least one axis but does not necessarily wrap around.) In Ref. NewmanZiff01,  we showed that many estimates of $`p_c`$ on the torus converge a factor $`L`$ faster than estimates for the open square—some converging as fast as $`L^{21/\nu }`$.

In conclusion, it is clear that the convergence of estimates for the critical occupation probability $`p_c`$ in percolation systems is highly dependent upon the nature of the estimate, as well as the shape and boundary conditions of the system, and that the shrewd use of this fact can allow one to make very accurate estimates of $`p_c`$ and scaling exponents.

## Appendix A Exact enumeration results

In Table 2 we give the exact expressions for the crossing probability function $`R_L(p)`$ for site percolation on a square lattice of size $`L\times L`$, for crossing from one given side of the square to the opposite side (such as left to right). The results for $`L=2`$ to 5 were given previously by Reynolds et al. ReynoldsStanleyKlein80 . Those for $`L=6`$ and $`L=7`$ were calculated previously for the work reported in Ref. Ziff92, , but reported in a different format (results were given for $`R_L^{}(p)`$ rather than $`R_L(p)`$ itself). From the results here, one can with reasonable ease calculate the various estimates of $`p_c`$ given in Table 1 using a symbolic manipulation program such as Maple or Mathematica. A file containing these polynomials in forms readable by such programs is available by email from the authors.

An alternative way to represent these results is as a series in $`p^nq^{Nn}`$ where $`q=1p`$ and $`N=L^2`$. The transformation can be achieved by substituting $`p1/(1+r)`$, multiplying by $`(1+r)^Np^N`$, expanding, and replacing $`rq/p`$ again. The results are given in Table 3.

This is also the form that Reynolds et al. ReynoldsStanleyKlein80  used in their series for $`R_2`$ through $`R_5`$. From the present point of view, these series are interesting because they are in precisely the form of Eq. (15), so that the coefficient $`c_{L,n}`$ of $`p^nq^{Nn}`$ above is related to the microcanonical crossing probability $`R_n`$ simply by

$$R_{L,n}=\frac{c_{L,n}}{\left(\genfrac{}{}{0pt}{}{N}{n}\right)}=\frac{c_{L,n}n!(Nn)!}{N!}.$$
(33)
That is, $`c_{L,n}`$ represents the number of configurations with $`n`$ occupied sites that satisfy the crossing criterion, out of a total of $`\left(\genfrac{}{}{0pt}{}{N}{n}\right)`$ possible configurations of the $`n`$ occupied sites among the $`N=L^2`$ sites of the lattice. For example, of the $`16!/(6!10!)=8008`$ possible configurations of 6 occupied sites on a $`4\times 4`$ lattice, exactly 390 are percolating by crossing in one direction (from the third term in $`R_4(p,q)`$), yielding a microcanonical probability $`R_{4,6}=390/8008=0.048701\mathrm{}`$ Likewise, for $`n=16`$ occupied sites on the $`4\times 4`$ system, there is exactly one percolating system out of one total system.

Thus, the polynomials given in Table 3 represent the microcanonical $`R_{L,n}`$ for $`L`$ up to 7.