text:
# Adsorption of polyampholytes on charged surfaces

## 1 Introduction

The study of the adsorption of charged polymers on charged surfaces is of great interest in both polymer physics and industrial practice. A better understanding of the adsorption phenomenon has already led to numerous applications related to wetting, lubrication, adhesion…Polyelectrolytes are for instance currently used to control the colloidal stability of dispersions or in waste water treatment. Relatively few experimental studies have however been performed on polyampholyte adsorption. Polyampholytes are polymers carrying both positive and negative charges along the same chain blaakmeer ; peffer ; kato ; leberre . On the theoretical side, Dobrynin, Rubinstein and Joanny have recently proposed a new adsorption mechanism induced by the polarization of the chains in the electric field of the surface dobrynin .

The present work aims at comparing theory and experiments by characterizing adsorbed polyampholyte layers on controlled charged surfaces with atomic force microscopy (AFM). Section 2 presents the experimental procedures: the fabrication of the charged surfaces by reaction of neutral and charged thiols on gold-coated substrates, the AFM experiments and the two characterization techniques of the surfaces, contact angle and streaming current measurements. Section 3 summarizes the theory of polyampholyte adsorption: the case of a single chain dobrynin  has been extended to the case of a solution obukhov  and the monomer concentration profile allows us to derive the theoretical loop size distribution in the adsorbed layer. In Section 4, the experimental loop size distribution is compared to theory and qualitative agreement is found. Contact angle and electrokinetics measurements are interpreted in the light of the AFM results.

## 2 Experimentals

### 2.1 Surface preparation

The charged surfaces have been prepared as described in bain1 . Substrates were cleaved mica for AFM experiments and contact angle measurements, or microscope glass slides for streaming current measurements. The gold-coated surfaces have been realized under high vacuum by deposition of 5 nm of chromium, which acts as a wetting agent, followed by 50 nm of gold. The surfaces are then left about 20 hours in an ethanol millimolar solution of thiols. We have used HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$ Na+, HS-(CH2)2-CH3 (both purchased from Aldrich) and HS-(CH2)2-(CF2)7-CF3 (kindly provided by Marie-Pierre Kraft from ICS-Strasbourg). The surface charge is changed by mixing charged and neutral thiols in varying proportions. Table 1 gives their respective lengths when the aliphatic (or fluorinated) chain is extended.

### 2.2 Polyampholyte adsorption

The polyampholyte has been synthesized by microemulsion-polymerization by S. Neyret microemul  from a neutral monomer, acrylamide, and two charged monomers, sodium 2-(acrylamido)-2-methylpropanesulfonate (AMPSNa) and \[2-(methacryloxy)ethyl\]trimethylammonium chloride (MADQUAT). The chemical formulae are shown on Figure 1.

The total number of monomers in the polymer is $`\mathrm{4\hspace{0.17em}10}^4`$; its composition (93.5/3.75/3.75 % in acrylamide/AMPSNa/MADQUAT) has been determined by elementary chemical analysis and shows that the polyampholyte is globally electrically neutral. The experimental determination of reactivity ratios has shown that polymerization in microemulsion leads to an almost random distribution of monomers along the chain microemul . This is due to the strong concentration of counterions in the microemulsion droplets (diameter $``$ 50-100 nm) that strongly screens the electrostatic interactions between charged monomers. However, because of the different reactivity ratios of the monomers, the fluctuations of the net charge of the polyampholytes (defined as the sum of all the charges along the chain), are larger than predicted for a perfectly random distribution diffusion . After the synthesis, the counterions are removed by dialysis against pure water. The polyampholyte solution has been prepared by dissolving 10 mg of polymer in 100 g of a NaCl aqueous solution (concentration 0.15 mol.l-1), the polymer being insoluble in pure water and the critical salt concentration to reach solubility, as determined in neyret , being 0.1 mol.l-1. The thiolated surfaces have been left in the polyampholyte solution for one hour, this was long enough for the polymer to adsorb. After that, the surfaces have been rinsed with milli-Q water before being used in the AFM experiments. While salt is needed in bulk solution to ensure that polyampholyte chains do not collapse, the electric field of a charged substrate should be sufficient to stretch adsorbed chains (cf Fig. 5) and avoid their collapse; loops of adsorbed polyampholyte can then develop in milli-Q (de-ionized) water. This has been addressed theoretically in netz  where it is suggested that the behaviour of an adsorbed polyampholyte is similar to the behaviour of an adsorbed soluble polymer. Throughout our study, we consider water as a $`\theta `$-solvent for the polyampholyte; water is indeed close to $`\theta `$-solvent for polyacrylamide (the second virial coefficient $`A_2`$ for a polyacrylamide with a molar mass equal to $`\mathrm{4.7\hspace{0.17em}10}^6`$ g.mol-1 is low: $`A_2\mathrm{0.64\hspace{0.17em}10}^4`$ mol.cm3.g-2 viriel ) and the used polyampholyte is mostly composed of acrylamide.

### 2.3 AFM experiments

A Digital Instruments Nanoscope III force microscope has been used, fitted with its standard fluid cell. For each experiment, a new cantilever was cleaned in a water plasma plasma . Cantilever spring constants (from Digital Instruments specifications) were 0.05 N.m-1 for experiments on HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$/ HS-(CH2)2-CH3 surfaces and 0.03 N.m-1 for experiments on HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$/ HS-(CH2)2-(CF2)7-CF3 surfaces. Milli-Q water (pH = 5.7) is injected in the liquid cell set-up and we have then waited half an hour to start the AFM force measurements. The tip-surface then underwent repeated approach-separation cycles at 1 Hz frequency; beyond this value, hydrodynamic forces appear. In de-ionized water (pH = 5.7), the total salt concentration (H3O+ and HCO$`{}_{}{}^{}{}_{3}{}^{}`$) is $`\mathrm{4\hspace{0.17em}10}^6`$ mol.l-1, which corresponds to a Debye length $`\kappa ^1150`$ nm. The thickness of the adsorbed layer predicted by theory being roughly $`\kappa ^1`$ in the pseudo-brush regime, a $`z`$-displacement of the piezoelectric system much larger than this length was chosen, namely 500 nm, which allows to determine precisely the force profile over a 400 nm range.

A peak sometimes appears on the force curves and is understood as the detachment of a loop that bridged tip and substrate senden ; silber (Figure 2). Some extension curves have been fitted using the Freely Jointed Chain fjc  and the Worm Like Chain wlc  models and the extension of the loop at the detachment was always found above 90 % of the contour length; the detachment distance of the loop can be considered as its contour length.

For each AFM experiment, more than two hundred jumps have been collected to obtain the loop size distribution. The data have been analyzed using the method introduced in senden  for the adsorption of polydimethylsiloxane on silica, and later used to investigate the adsorption layers of polyacrylamide derivatives silber . The loop size distribution $`S(n)`$ (defined as the probability for a loop to have $`n`$ monomers) is deduced from the distribution $`p(D)`$ of detachment distances $`D`$ via the relation $`S(n)dn=p(D)dD`$ and the hypothesis that $`Dna`$ where $`a`$ is the Kuhn length. Assuming that the experimental situation corresponds to the fence regime defined in Section 3 and using Eq. (13), the distribution $`p(D)`$ is fitted by:

$$p(D)=AD^{1/2}(D+h)^1$$
(1)
$`A`$ is a numerical prefactor and the number $`g`$ of monomers in a chain section with a size of the order of the layer of thickness $`\lambda `$ (the Gouy-Chapman length, see Section 3) is given by $`h/a`$. Values of $`D`$ lower than 20 nm have not been taken into account in the fit because of a lower resolution near the surface and the presence of a primary adhesion peak (due to van der Waals interactions) on some force curves.

### 2.4 Contact angle measurements

Contact angles of a water sessile drop have been measured using an optical microscope. The surfaces were tilted until the drop moved, defining both advancing $`\theta _a`$ and receding $`\theta _r`$ angles. Because of the Young-Dupré relationship $`\gamma _{LV}\mathrm{cos}\theta =\gamma _{SV}\gamma _{SL}`$, where $`\gamma _{LV}`$, $`\gamma _{SV}`$ and $`\gamma _{SL}`$ are respectively the interfacial free energies of the liquid-vapour, solid-vapour and solid-liquid interfaces, $`\mathrm{cos}\theta `$ was plotted against the ratio of charged thiols to (charged + uncharged) thiols. Indeed, for a given liquid/vapour couple (here water/air), $`\mathrm{cos}\theta `$ should approximately vary linearly with the solid surface composition if there is neither segregation of the thiols at the surfaces nor special conformational effects of the thiols. The contact angle hysteresis $`(\mathrm{cos}\theta _r\mathrm{cos}\theta _a)`$ also gives some insight on the surface heterogeneities hyster1 ; hyster2 ; hyster3 .

### 2.5 Electrokinetic measurements

An original apparatus to measure the streaming current of planar substrates has been designed and is schematically represented on Figure 3. Figure 4 shows a detailed view of the measurement cell.

The measurement cell is built from two PMMA plastic blocks that comprise an accurately cut depression to fit in two test glass slides of dimensions 7.5 x 2.5 x 0.1 cm3. Each slide is firmly held in position by two Viton rubber stops. A small reservoir was machined at each extremity of the base block. The two blocks are separated by a Teflon (PTFE) gasket of 0.02 cm thickness, aligned with stainless steel pegs. A capillary of dimension 7.5 x 2.2 x 0.02 cm3 is then formed between the two slides by clamping together the two blocks with 8 screws. The principle of the measurement is the following: upstream, an electrolyte solution ($`10^3`$ mol.l-1 KCl at pH = 5.7) falls into a column so that the level of the liquid remains constant, which ensures a constant pressure difference $`\mathrm{\Delta }p`$ between the two extremities of the measurement cell; $`\mathrm{\Delta }p`$ is measured via two Kobold pressure transducers (0-$`10^4`$ Pa) linked to tightly closed reservoirs. The measured pressure difference accuracy is 5 Pa. The streaming current is measured with a Keithley 617 electrometer (impedance $`\mathrm{2\hspace{0.17em}10}^{14}`$ M$`\mathrm{\Omega }`$) used as a nano-amperemeter connected to two silver electrodes placed in each reservoir of the measurement cell. The electrodes were covered with silver chloride by electrolysis to prevent their polarization.

The principle of the $`\zeta `$ potential determination has been discussed by Hurd and Hackerman hurd  in their study of metals (Au, Ag and Pt) in aqueous solutions. In order to measure the streaming current, we have used a shunt resistance of $`10^3`$ $`\mathrm{\Omega }`$; the measured impedance of the capillary tube is 1 M$`\mathrm{\Omega }`$ which ensured that the measured current did not depend on the value of the shunt resistance hurd .

In the case of a laminar and established flow across a charged capillary of length $`L`$, width $`l`$ and height $`h`$, the $`\zeta `$ potential can be calculated from the pressure $`p`$ dependence on the streaming potential $`U_s`$ and the streaming current $`I_s`$ (Smoluchowski equations hunter ):

$$\zeta (U_s)=\frac{\eta \left[\lambda _l+{\displaystyle \frac{2\lambda _s}{h}}\right]}{\epsilon }\frac{dU_s}{d(\mathrm{\Delta }p)}$$
(2)
$$\zeta (I_s)=\frac{\eta }{\epsilon }\frac{L}{lh}\frac{dI_s}{d(\mathrm{\Delta }p)}$$
(3)
where $`\zeta (U_s)`$ is the $`\zeta `$ potential calculated from the streaming potential, $`\zeta (I_s)`$ the $`\zeta `$ potential obtained from the streaming current, $`\eta `$ the viscosity, $`\epsilon `$ the permittivity of the fluid, $`\lambda _l`$ the liquid (electrolyte) conductivity and $`\lambda _s`$ the surface conductivity. If the influence of the interfacial conductivity is not considered, only an apparent $`\zeta `$ potential is determined from the streaming potential which is lower than the actual one. The $`\zeta `$ potential determined from the streaming current, however, is independent of the surface conductivity and is a function only of the channel geometry for a given electrolyte. The thiolated gold surfaces are a priori conducting and we have then chosen to deduce $`\zeta `$ from streaming current measurements.

## 3 Theory

We summarize in this section the theory of polyampholyte adsorption in a $`\theta `$-solvent dobrynin ; obukhov  that we need to derive the loop size distributions. We also give numerical estimates of the various threshold values using the characteristics of the polyampholyte chains used in our experimental study and described in Section 2 ($`N`$, number of monomers $`\mathrm{4\hspace{0.17em}10}^4`$; $`f`$, fraction of charged monomers $`0.075`$; $`a`$, monomer size $`0.3`$ nm; $`R_0`$, ideal (Gaussian) chain size $`60`$ nm).

### 3.1 Single chain adsorption dobrynin

In the Poisson-Boltzmann approximation, a charged surface (of surface charge density $`\sigma _0`$), immersed in a solvent of dielectric permittivity $`\epsilon `$, in absence of added salt, generates an electric field $`E(z)`$ at a distance $`z`$ from the surface:

$$E=\frac{e\sigma _0}{\epsilon (1+z/\lambda )}$$
(4)
where $`\lambda `$ is the so-called Gouy-Chapman length which gives the thickness of the counterion layer close to the surface israelachvili :

$$\lambda =\frac{1}{2\pi \sigma _0l_B}$$
(5)
with $`l_B`$ the Bjerrum length $`l_B=e^2/(4\pi \epsilon k_T)`$ ($`e`$ is the elementary charge and $`k_BT`$ the thermal energy, $`l_B7`$ Å in pure water).

Because of the statistical distribution of the charges along the chain, the two halves of a neutral chain carry opposite charges of order $`(fN)^{1/2}e`$. The electric field can polarize a neutral polyampholyte chain which behaves as a dipole (of size $`R_z`$ and charge $`e(fN)^{1/2}`$) and is then attracted by the charged surface. In a $`\theta `$-solvent of the polymer backbone, the chain is stretched when $`R>N^{1/2}a`$; this corresponds to a threshold field $`E_11/(N^{1/2}a)(k_BT/e)(fN)^{1/2}`$.

Three regimes of adsorption have been predicted depending on the charge density $`\sigma _0`$ (Figure 5).

i) pole regime ($`\lambda aN^{1/2}`$). A chain is stretched when $`E>E_1`$; the corresponding surface charge density and Gouy-Chapman length are respectively $`\sigma _1(fN)^{1/2}(l_BR_0)^1`$ and $`\lambda _1R_0(fN)^{1/2}.`$ For the polyampholyte studied in the experiment, we get $`\sigma _1\mathrm{6.9\hspace{0.17em}10}^{13}`$ charges/m2 and $`\lambda _13290`$ nm. Beyond this threshold, the chain size $`R_z`$ in the direction normal to the surface grows linearly with charge density as $`R_0\sigma _0/\sigma _1.`$

ii) fence regime. The Gouy-Chapman length $`\lambda `$ becomes comparable with $`R_z`$ when $`\sigma _0\sigma _2(fN)^{1/4}(l_BR_0)^1`$. Because of the screening of the electric field by the counterions, the chain then remains confined within a slice of height $`\lambda `$ (for our system, $`\sigma _2\mathrm{5.2\hspace{0.17em}10}^{14}`$ charges/m2 associated with a Gouy-Chapman length $`\lambda _2440`$ nm). An adsorbed chain can be divided into subunits (blobs) of size $`\lambda `$. The number $`g`$ of monomers per blob is obtained by considering each subunit as an independent chain:

$$gf^{1/3}(l_Ba\sigma _0)^{4/3}$$
(6)
Each blob carries an induced dipole moment $`e\lambda f^{1/2}g^{1/2}`$.

iii) pancake regime. The polyampholyte chain is strongly bound to the surface when the adsorption energy of each blob of size $`\lambda `$ is larger than the thermal energy $`k_BT`$. This happens for a charge density $`\sigma _0`$ larger than $`\sigma _3f^{1/2}/(al_B)`$. The condition $`\sigma _0>\sigma _3`$ then defines the pancake regime: the Gouy-Chapman length $`\lambda `$ is smaller than the mean-square distance between charged monomers $`af^{1/2}=\lambda _3`$ ($`\sigma _3\mathrm{2.3\hspace{0.17em}10}^{17}`$ charges/m2 and $`\lambda _31`$ nm in our experimental situation). Monomers with charge opposite to the surface charge are in contact with the surface while the other monomers form loops dangling in solution at distances $`z>\lambda `$. The average size of these loops can be estimated from the balance between the electrostatic force acting on the monomers $`e^2\sigma _0\lambda /(\epsilon z)k_BT/z`$ and the elastic force $`k_BTzf/a^2`$ needed to stretch a chain with $`f^1`$ monomers; the thickness of the adsorbed layer is then approximately $`af^{1/2}`$. It is interesting to note that this thickness does not depend on the surface charge $`\sigma _0`$.

### 3.2 Multichain adsorption obukhov

The previous model has been extended to the case of multichain adsorption by Dobrynin and co-workers obukhov .

i) multilayers of stretched chains ($`\sigma _1<\sigma _0<\sigma _2`$, Figure 6.a).

The predicted monomer concentration profile is:

$$c(z)c_0^{}\frac{\lambda _1}{\lambda +z}$$
(7)
where $`c_0^{}=a^3N^{1/2}`$ is the overlap concentration. At distances $`z>\lambda `$, a hyperbolic density profile $`c(z)c_0^{}\lambda _1/z`$ is expected. Near the surface, the polymer density saturates at $`c(0)c_0^{}\lambda _1/\lambda `$. The crossover between semi-dilute and dilute regimes of the adsorbed chains occurs at distance $`z\lambda _1`$ where the attractive interaction between each chain and the surface compares to $`k_BT`$: $`\lambda _1`$ can be considered as the thickness of the adsorbed layer.

ii) self-similar stretched pseudo-brush ($`\sigma _2<\sigma _0<\sigma _3`$, Figure 6.b).

The size of the chains in the first adsorbed layer saturates at $`\lambda _2R_0(fN)^{1/4}`$. The adsorbed layer at distances $`z<\lambda _2`$ can be viewed as a self-similar pseudo-brush of stretched polydisperse loops.

The stretching of a blob of size $`z`$ containing $`g(z)`$ monomers can be estimated by balancing the elastic energy $`k_BTz^2/(a^2g(z))`$ with the polarization energy $`e\sqrt{fg(z)}zE(z)`$. For $`z>\lambda `$ the electric field is $`E(z)e/(\epsilon l_Bz)`$. The equilibrium density profile is given by:

$$c(z)a^3\left(\frac{af}{z+\lambda }\right)^{1/3}$$
(8)
For $`z>\lambda _2`$, the concentration profile is still given by equation (7). The total thickness of the adsorbed layer is $`\lambda _1`$.

iii) $`\sigma _0>\sigma _3`$. The size of the polyampholyte molecules in the first layer is $`\lambda _2`$, but these chains are self-similarly stretched at all length scales between $`\lambda _3`$ and $`\lambda _2`$. In the layer of thickness $`\lambda _3`$ near the wall the monomer concentration is constant and proportional to $`a^3f^{1/2}`$. At distances between $`\lambda _2`$ and $`\lambda _1`$ one finds a multibrush of stretched chains.

### 3.3 Loop size distribution inside the pseudo-brush for $`\sigma _2<\sigma _0<\sigma _3`$

We derive in this section the loop size distribution in the regime of interest for our experiment. The same treatment can be used to obtain the loop size distribution in the other regimes.

We consider a layer of thickness $`dz`$ located at a distance $`z`$ from the charged surface; the number of monomers (per unit area) in this layer is degennes :

$$c(z)dz=p(n)dn$$
(9)
where $`c(z)`$ is the monomer concentration and $`p(n)`$ the probability (per unit area) that a monomer belongs to a loop comprising more than $`n`$ monomers. The loop size distribution $`S(n)`$, probability to find a loop with n monomers, is then derived from $`p(n)`$ by:

$$S(n)=\frac{p}{n}$$
(10)
The relation between $`z`$ and $`n`$ is found by balancing the elastic energy of a strand of size $`z`$ containing $`n`$ monomers and its polarization energy in the electric field $`E(z)`$:

$$\frac{k_BT}{na^2}z^2\frac{e\sigma _0}{\epsilon }\frac{\lambda }{z+\lambda }e(fn)^{1/2}z$$
(11)
leading to:

$$z(z+\lambda )n^{3/2}$$
(12)
Using equation (8) for $`c(z)`$, we obtain $`S(n)\lambda ^{4/3}n^{1/2}`$ for $`z<\lambda `$, and $`S(n)n^{3/2}`$ for $`\lambda <z<\lambda _2`$. These two asymptotic expressions of $`S(n)`$ can be recast into a single formula:

$$S(n)n^{1/2}(g+n)^1$$
(13)
where $`g`$, defined by equation (6), is the mean number of monomers in a strand of size $`\lambda `$. The polyampholytes located at distance $`z>\lambda _2`$ are not adsorbed directly on the surface and should not be stretched in the AFM experiments; they thus do not participate to the loop size distribution. This loop size distribution 13 is identical to the one derived by Dobrynin et al dob  when $`n>g`$; but for $`n<g`$, we find a power-law profile instead of a constant loop size, although we have to assume that all chains fold back for $`z=\lambda `$.

## 4 Results and discussion

Images obtained by AFM on the HS-(CH2)2-(CF2)9-CF3 monolayer showed a remarkably ordered 2D hexagonal lattice with a lattice constant equal to 0.59 nm on Au(111)/mica tamada ; the surface density is then $`\mathrm{3.3\hspace{0.17em}10}^{18}`$ thiols/m2. Keeping the same value for our systems and making the assumption that the ratio of charged thiols to (charged + neutral) thiols on the gold surface remains equal to the ratio of charged thiols to (charged + neutral) thiols in the bulk solution used to process the gold surface (a very strong assumption indeed), we can estimate the surface charge density $`\sigma _0`$ for the various substrates. Using the estimated values of $`\sigma _1`$, $`\sigma _2`$ and $`\sigma _3`$ of Section 3, it appears that our experiments are performed in the fence regime. This is justified later on Figure 11 and allows us to process our data according to the theoretical analysis exposed in 3.3. The last assumption supposes that there is no differential adsorption of the two thiols on gold.

### 4.1 Loop size distribution on HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$/HS-(CH2)2-CH3

The loop size distribution is shown on Figure 7.

We find from the numerical fit that $`g0`$ which means that no monomers are adsorbed in the layer of thickness $`\lambda `$. The loop size distribution $`S(n)`$ scales as n-3/2 as in the studies of neutral polymers reported in senden ; silber . In fact, we have checked that the polyacrylamide homopolymer adsorbs on a surface covered by HS-(CH2)2-CH3; short-chain thiols form less dense monolayers than long-chain thiols bain1  and are not efficient enough to prevent the polyampholyte adsorption; short-chain neutral thiols are thus not good candidates to investigate electrostatic effects.

### 4.2 Loop size distribution on HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$/HS-(CH2)2-(CF2)7-CF3

Giving up with HS-(CH2)2-CH3, we turned to HS-(CH2)2-(CF2)7-CF3 because of its longer chains and enhanced hydrophobicity. We have first checked that polyacrylamide does not adsorb on the fluorinated thiols: no detachment peak has been observed on the force curves. We cannot of course rule out from this observation that the polyampholyte cannot adsorb by non-electrotatic interactions between the fluorinated thiols and the chain skeleton; this possibility has not been investigated. On the other hand, polyacrylamide adsorbs on a surface covered with HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$, certainly via van der Waals interactions with the underlying gold surface, and consequently the polyampholyte also adsorbs. In order to favour the electrostatic interactions between the charged monomers and the surface against the hydrogen bonds between acrylamide and the charged thiols, weakly charged surfaces have been prepared with millimolar solutions of thiols mixtures HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$ Na+/HS-(CH2)2-(CF2)7-CF3 with molar ratios 0.1, 0.3, 1 and 3 %. For the two extreme cases, the loop size distribution and the fit given by equation (1) are shown on Figures 8 and 9.

For each type of surface (0.1, 0.3, 1 and 3 %), experiments have been repeated with four different substrates and the results of the fit (using Eq. 1) are given in Table 2. The different results obtained show that there is a charge effect. We can deduce that the salt initially added to dissolve the polyampholyte chains has been effectively removed from the adsorbed layer after rinsing, otherwise the field at the surface would be screened after less than 1 nm (for an electrolyte concentration of 0.15 mol.l-1) and then the loop size distributions for different surfaces would be identical. Moreover, the preferential adsorption of one thiol in comparison with the other is limited, otherwise the surfaces would be covered with the same type of thiol and there would not be any charge effect. This observation does not however allow to conclude that the rate of charged thiols at the surface is the same as the one in volume.

The results for the mean number $`<g>`$ of monomers in a loop are also plotted on Figure 10 taking the standard deviation as experimental bars.

Figure 11 shows that the AFM experiments are in the fence regime; as predicted by theory, the number $`g`$ of monomers in a loop decreases with the surface charge. $`g`$ scales like the ratio of charged/neutral thiols to the -0.86 power; however, it seems difficult to compare it to the exponent -1.33 given by Eq. 6, because the theoretical exponent links $`g`$ to the surface charge and not the bulk fraction of charged thiol.

### 4.3 Contact angle and electrokinetics measurements

Contact angles and $`\zeta `$ potential for the charged surfaces made with HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$/HS-(CH2)2-(CF2)7-CF3 are presented on Figure 12.

The values of $`\zeta `$ obtained for the monolayer of fluorinated thiols (- 37 and - 39 mV) are close to the value measured by Chibowski and Waksmunki for polytetrafluoroethylene (CF2)n teflon . In bidistilled water, they indeed found a potential equal to - 46.6 mV. The $`\zeta `$ potential of the monolayer of charged thiols, - 55 mV, can be compared to that of a surface covered by mercaptopropionic acid (HS-(CH2)2-COOH). By atomic force measurements in a $`10^3`$ mol.l-1 KCl aqueous solution, Hu and Brad bard  determined the potential $`\psi _d`$ above the counterion layer of this surface; $`\psi _d`$ is indeed close to $`\zeta `$ hunter . At pH $`>10`$, the acid is completely dissociated, and the potential $`\psi _d`$ was - 62 mV, which is compatible with the streaming current measurement.

Figure 12 shows that $`\mathrm{cos}\theta `$ and $`\zeta `$ are relatively well correlated: their values are hardly modified for a bulk fraction of charged thiol lower than 70-75 %, and then strongly change. The contact angle and the $`\zeta `$ potential are thus linked to the same phenomenon. The plateau can be explained by i) a preferential adsorption of the fluorinated thiol, ii) a phase separation between the two types of thiols, or iii) the conformations of the thiol alkyl chains. The first hypothesis is refuted by the AFM experiments, which show a difference between the surfaces in the region of the plateau (0.1 - 3 %). We have not performed AFM imaging to investigate the second hypothesis and the possible presence of islands resulting from phase separation. Anyway the contact angle hysteresis $`\mathrm{cos}\theta _r\mathrm{cos}\theta _a`$ is limited and identical to the one of a pure monolayer of fluorinated thiols, which is a very good indication of an homogeneous surface; the $`\zeta `$ potential would also be extremely sensitive to the presence of charged islands. The third hypothesis seems to be the most likely. As shown by Figure 13, the HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$ thiols are much smaller than the fluorinated ones, and the charges that they carry are hidden by the long neutral chains, which decreases their influence by a screening effect: the local reduced dielectric constant should also decrease the dissociation constant of the sulfonic groups of the smaller thiols, with the effect of a smaller effective smaller charge. This remains very schematic, since thiolated chains can be tilted on the surface and the long chains are unlikely to be in an all-trans conformation. However, AFM images of a monolayer of the fluorinated thiol have revealed extended conformations with a tilt angle of 22 naud .

Figure 14 represents data obtained by A.-L. Bernard, who studied the surfaces HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$/HS-(CH2)2-CH3 annelaure . Contrary to the former situation, the variation of $`\mathrm{cos}\theta `$ with respect to the percentage of charged thiols in the bulk is nearly linear, except for highly charged surfaces. The two thiolated molecules have the same length and there is no screening effect.

### 4.4 Interpretation

A striking feature of Figure 12 is the non-zero value of the $`\zeta `$ potential of the monolayer of fluorinated thiols, which are neutral molecules. Similar streaming current measurements on a monolayer of octadecanethiol zetathiol  have showed that the $`\zeta `$ potential varies with pH and salt concentration. The proposed explanation is a preferential adsorption of the hydroxide ions on the surface. In the case of the fluorinated monolayer, a preferential adsorption of anions certainly takes place, but it is difficult to decide whether the adsorbing ions are OH- or Cl-. The phenomenon is sketched on Figure 15. The $`\zeta `$ potential is approximately equal to the potential $`\psi _d`$ at the outer Helmoltz plane (OHP) and thus non-negligible due to the adsorption of anions, whereas the surface potential $`\psi _0`$ remains weak hunter .

The surface charge $`\sigma _d`$ of the OHP plane close to the monolayer of charged thiols can be obtained with the hypothesis $`\zeta \psi _d`$ and the relation given by the Gouy-Chapman-Stern-Graham theory: $`\sigma _d=4ce\kappa ^1\mathrm{sinh}(e\psi _d/2k_BT)`$, where $`c`$ is the salt concentration, and $`\kappa ^1`$ the Debye length hunter . In our streaming current measurements on a monolayer of charged thiols, c = $`10^3`$ mol.l-1, $`\kappa ^1`$ 10 nm, $`\psi _d`$ = - 55 mV. The obtained value $`\sigma _d`$ $`\mathrm{4.9\hspace{0.17em}10}^3`$ C.m-2 corresponds to a charge density equal to $`\mathrm{3.1\hspace{0.17em}10}^{16}`$ negative charges / m2, much lower than the surface density expected for a monolayer of charged thiols: low energy helium diffraction on a HS-(CH2)9-CH3 monolayer helium  and surface acoustic wave studies on a HS-(CH2)6-CH3 monolayer acoustic  both established that the surface density of the molecules of the thiolated monolayer is around $`\mathrm{5\hspace{0.17em}10}^{18}`$ thiols/m2. If we suppose that the density is the same for the charged thiol monolayer, the large difference between this value and the former charge density implies a strong screening by the counterions (cations) in solution, represented on Figure 16. A similar effect was also inferred by Hu and Hard after AFM experiments on thiol monolayers hu . We would like to stress here that it is anyway difficult to determine the actual charge that rules the electric potential of a charged surface. The determination of the charge is in fact model-dependent; the actual charge density is likely to be between the values estimated from direct (crystallographic) measurements (such as AFM tamada ) and electrokinetic measurements.

Concerning the surfaces used in AFM experiments (0.1, 0.3, 1 and 3 % of charged thiols in the bulk), contact angle and electrokinetics measurements do not show any difference with respect to a monolayer of 100 % neutral fluorinated thiols (see Figure 12): $`\theta 115^{}`$ and $`\zeta 38`$ mV. The study of the detachment distances in AFM experiments have however revealed a charge effect in qualitative agreement with the theory. These apparently contradictory observations can be explained by making the following assumption: the few short-chain charged thiols are hidden by the long fluorinated ones. Contact angle and electrokinetics measurements give then the feeling that these surfaces are covered by a monolayer of neutral thiols. These measurements are indeed realized at a mesoscopic scale, whereas the AFM ones are realized at a microscopic scale. Even if the contact angle and the zeta potential see only the neutral fluorinated thiols, polyampholytes have the ability to find their way between them to adsorb on the charged thiols (Figure 17).

## 5 Conclusion

The adsorption of a polyampholyte has been studied on various HS-(CH2)2-SO$`{}_{}{}^{}{}_{3}{}^{}`$/HS-(CH2)2-(CF2)7-CF3 monolayers by AFM. The statistical analysis of the loop size distribution in the adsorbed layer reveals that the number of monomers in a loop decreases, in qualitative agreement with theoretical predictions dobrynin ; obukhov . In order to measure the charge density $`\sigma _0`$ of the surfaces, we have performed streaming current measurements, which gave us the $`\zeta `$ potential of the surfaces. Because of a phenomenon of preferential adsorption of some ions, the value of $`\sigma _0`$ could not be deduced from $`\zeta `$. For instance, the existence of a non-zero $`\zeta `$ potential on the monolayer of fluorinated thiols is difficult to interpret; an explanation could be that the bond between the fluorocarboned part of the thiol and its hydrocarboned part creates a strong dipole dipole ; this could attract anions of the electrolyte by induced-type interactions. The dipoles oriented in the bonds -CH2-CF2\- moreover constitute a layer of dipoles which could contribute to this attraction.

The contact angle and $`\zeta `$ potential measurements are well correlated, and give the same value for the monolayer of fluorinated thiols than for weakly charged surfaces. AFM experiments on these surfaces however show different structures of the polyampholyte adsorbed layer. This apparent paradox is explained by the following interpretation: the short-chain charged thiols are hidden by the long-chain fluorinated ones, but the polyampholyte has the ability to find its way between them to adsorb. This hypothesis could be checked by infrared or Raman surface spectrophotometry and, if confirmed, could lead to interesting applications, such as the recognition by charged polymers of charged groups which are undetectable at a mesoscopic scale by a simple wetting or electrokinetics study.

## 6 Acknowledgments

The authors are thankful to Drs. Françoise Candau, Franck Clément, Georges Debrégeas, Albert Johner, Marie-Pierre Kraft, Hans Lyklema, Philippe Mesini, Joseph Selb and Tim Senden for fruitful discussions.


format: text

text:
# The 8Li Calibration Source For The Sudbury Neutrino Observatory

## 1 Introduction

The Sudbury Neutrino Observatory (SNO) is a 1 kilotonne mass heavy water Cherenkov detector designed for detection of neutrinos from the 8B reaction in the Sun. Neutrinos interact in D2O by elastic scattering on electrons (ES), and the charged current (CC) and neutral current (NC) interactions on deuterons. The ES and CC interactions both result in Cherenkov light from an energetic electron, and the NC interaction results in the disintegration of the deuteron which is observed by $`\gamma `$-rays from the capture of the free neutron. Cherenkov light from these interactions is detected by a spherical array of $``$ 9500 photomultiplier tubes surrounding the D2O volume, which is contained in a 12 meter diameter acrylic vessel. The detector is immersed in $``$7 kilotonnes of light water to provide support and shielding. The detector is located in a cavity in INCO’s Creighton mine near Sudbury, Canada. A general description of the SNO detector and its subsystems can be found in Reference .

Cherenkov light measurements in SNO are based upon coincident phototube signals registered within an event (i.e. within 400 ns). Energy information is largely associated with the number of tubes hit (Nhit), and the event position can be reconstructed using the relative phototube timings within the event. The hit pattern can be used to determine direction information and discriminate against instrumental backgrounds. Use of these signals requires an understanding of the detector response both through the use of Monte Carlo techniques, and through analysis of calibration sources (also discussed in Reference ).

The SNO detector is calibrated using a variety of sources, including an isotropic light source, $`\gamma `$-rays, neutrons, and $`\beta `$-particles. This work describes a source which provides a high-endpoint energy spectrum of $`\beta `$-particles created by the decay of 8Li. The source can be used to assist in energy calibrations, position reconstruction studies, and evaluation of instrumental background discrimination.

The 8Li isotope is of interest for several reasons. The isotope decays with a Q-value of 16.003 MeV and the main beta-decay branch is an allowed decay with a central end-point energy of 12.96 MeV . The energy spectrum produced is similar to that of the 8B charged-current neutrino signal,  both in endpoint and shape, as 8Li $`\beta ^{}`$-decay is the isospin mirror to the 8B $`\beta ^+`$-decay (see Figure 1). Both isotopes decay to the same broad excited state of 8Be which in turn promptly decays to $`\alpha `$-particles that can be used to identify (“tag”) an event in the 8Li source.

The 8Li calibration source  consists of several subsystems, illustrated in Figure 2 and described in this paper. Radioisotope production and transport are described in sections 2 and 3. (Proof-of-principle and design tests of these mechanisms have been reported seperately ). The design of the decay chamber is described in section 4. The decay chamber is positioned using the SNO manipulator system, described in Reference . The performance of the source and its application are discussed in section 5.

## 2 8Li Production

The 8Li isotope is produced by the 11B(n,$`\alpha `$) reaction (Q=-6.63 MeV, $`\sigma `$=30.5 mb ) using a deuterium-tritium (DT) generator to provide approximately 108 14-MeV neutrons per second. This DT generator (MF Physics model A-320) and its environment are described more thoroughly in Reference . The primary production target consists of $`>99\%`$ enriched 11B boron powder suspended in a polyethylene matrix. The target lines the inner wall of an annular target chamber that surrounds the DT generator (see Figure 3). The outer wall of this chamber is lined with a polyethylene collar whose inner surface is coated with boron nitride.

The 8Li atoms recoil from the target into a laminar gas stream moving through the chamber. Helium is used as the transport gas, with sodium chloride aerosol added to provide an affinity for the charged lithium atoms. The aerosol is introduced to the gas stream using a tube oven containing a quartz tube with NaCl heated to approximately 610C, just below its melting point.

The chamber dimensions have been optimized to allow the recoiling 8Li atoms just enough room to stop in the gas before hitting a wall, but minimize the gas volume in the chamber so that the residence time of the gas is short. Because 8Li recoils in the direction of the incident neutrons, it is primarily the inner target which is responsible for 8Li production.

Although the total rate of isotope production in the target can be readily calculated, the recoil and uptake by the salt aerosol are not accurately predicted. Direct measurement of the collection rate is complicated because corrections based on the gas transport efficiency (see following section) are needed. Because of the many assumptions made in the gas transport model, it is difficult to separate 8Li collection from transfer efficiency when measuring yields. If the gas model described below is assumed correct then the 8Li recoil and capture rates (i.e. collection rate) can be estimated by correcting yield data with measurements from different transfer configurations. The estimated rate at which 8Li atoms are collected is (3$`\pm `$1) $`\times 10^7`$ 8Li atoms per 14-MeV neutron.

## 3 Gas Transport System

The isotope is transported from the shielded neutron generator pit to the detector by two sets of transfer lines. One set, a conduit approximately 40 meters long, connects the production chamber in the pit to the source manipulator system on a deck above the SNO detector. The second set is an umbilical approximately 26 meters long that runs through an umbilical retrieval mechanism, and down into a decay chamber immersed in the D2O. This umbilical length is sufficient to deploy the chamber in most of the top half of the D2O volume.

The construction of the umbilical is similar to that used for the SNO 16N source  but with larger overall dimensions. The design consists of a 3.38 mm inner diameter PTFE capillary used to flow gas into the decay chamber, and an 8.0 mm inner diameter polyethylene tube that is used to exhaust the gas from the decay chamber. The inlet capillary is inside the outlet tube. A signal line and low-voltage DC power lines are wrapped in a helix around the outer tube, and the whole assembly is covered with a clean, flexible silicone sheath acting as a water seal. Silicone potting gel is used to seal the space between sheath and outlet tube. This design allows robust manipulation of the umbilical and protects the radioactively pure D2O against contamination.

Because of the low production rate and short half-life of the 8Li isotope, the characteristics of the gas system (i.e. gas pressures, flow rates, capillary dimensions) were designed for optimization of the 8Li yield. The gas flow was designed to operate at the limit of laminar flow through the conduit and umbilical. Faster delivery of the radioisotope is possible through the use of near-supersonic gas flow speeds, but is deemed impractical in this case due to considerations of safety, complications such as salt buildup inside the decay chamber, and the difficulties of high pressure engineering.

The yield of 8Li may be taken as the product of the rate of radioisotope production ($`R_p`$), the efficiency with which the isotope can be extracted from the production chamber ($`ϵ_p`$), the efficiency for transport ($`ϵ_l`$), and the fraction of atoms that decay while inside the decay chamber ($`ϵ_d`$). These properties are estimated using gas transport dynamics as described in References .

The production chamber was designed for approximately 200-250 kPa absolute gas pressure to optimally range out 8Li recoils. Although the chamber was designed to promote laminar flow of gas, calculations were done for the pessimistic case of full mixing. The rate of change of the number of atoms of 8Li in the production chamber is described by Equation (1) for the case of full mixing,

$$\frac{dN}{dt}=R_p\lambda _{8Li}NE$$
(1)
where $`N`$ is the number of 8Li atoms in the chamber, $`R_p`$ is the rate of isotope pickup on the salt in the chamber, $`\lambda _{8Li}`$ is the decay rate of 8Li (ln2/t1/2), and E is the rate of 8Li atoms flushed from the chamber, such that $`E=\frac{N}{V_p}\frac{dV}{dt}`$, with $`V_p`$ as the volume of the production chamber.

We define the mass flow rate as $`Q=P\frac{dV}{dt}`$. The mean time to completely flush the production chamber, $`\tau _p`$, is determined by this flow rate.

$$\tau _p=\frac{1}{\lambda _p}=V_p\left(\frac{dV}{dt}\right)^1=\frac{PV_p}{Q}$$
(2)
The steady state solution of Equation (1) allows us to define the efficiency as the fraction of produced 8Li which escapes the target chamber, $`ϵ_p`$.

$$ϵ_p=\frac{E}{R_p}=\frac{\lambda _p}{\lambda _p+\lambda _{8Li}}=\frac{1}{1+\frac{\tau _p}{\tau _{8Li}}}$$
(3)
The mean time for the gas to move from the production chamber through the transfer line to the decay chamber, $`\tau _l`$ is:

$$\tau _l=_0^l\frac{P(x)}{Q}A𝑑x$$
(4)
where $`P(x)`$ is the pressure at point $`x`$ along the line, $`A`$ is the cross-sectional area of the transfer line, and $`l`$ is the total distance. For laminar flow this is calculable using Poisseuille’s law for mass flow, Equation (5). In this equation, $`\eta `$ is the gas viscosity, $`d`$ is the diameter of the transfer line, and $`P_p`$ and $`P_d`$ are the pressures of the production and decay chambers respectively.

$$Q=\frac{\pi d^4\left(P_p^2P_d^2\right)}{256\eta l}$$
(5)
For given pressures $`P_p`$ and $`P_d`$, Equation (4) can be computed to give the mean transfer time, Equation (6).

$$\tau _l=\frac{2Al}{3Q}\times \frac{P_p^3P_d^3}{P_p^2P_d^2}$$
(6)
The fraction of 8Li atoms that reach the far end defines the transfer line efficiency, Equation (7).

$$ϵ_l=\mathrm{exp}\left(\frac{\tau _l}{\tau _{8Li}}\right)$$
(7)
The decay chamber efficiency is the complement of the production chamber efficiency. Maximal gas mixing is assumed (now the optimal case) so that the efficiency is given by Equation (8).

$$ϵ_d=1\frac{\lambda _d}{\lambda _d+\lambda _{8Li}}=\frac{1}{1+\frac{\tau _{8Li}}{\tau _d}}$$
(8)
where $`\tau _d=\frac{1}{\lambda _d}=P_dV_d/Q`$, given the pressure and volume of the decay chamber as $`P_d`$ and $`V_d`$ respectively.

The net yield of the system is computed as the product of the above efficiencies.

$$\mathrm{Yield}=R_p\times ϵ_p\times ϵ_l\times ϵ_d$$
(9)
The system design is constrained by both practical considerations and by the desire to maximize this yield. To begin with, the volume of the decay chamber was chosen in order to get the best signal (see Section 4). The decay chamber pressure was determined by the conductance of the umbilical outlet tube, which was exhausted to laboratory atmospheric pressure, 125 kPa (absolute) at 2000 m below surface.

Arbitrarily high gas flow rates could not be used, because the resulting turbulent flow would not be efficient in transporting the radioisotope. Turbulent currents could increase the time it takes for a given particle to transit the transfer line, and could also deposit salt aerosol on the walls of the transfer line, both of which impair efficiency. Laminar flow occurs only when the flow rate is such that the Reynold’s number is less than a critical value of approximately 1000 . The Reynold’s number, $`R_e`$ is a dimensionless quantity defined as:

$$R_e=\frac{4m}{\pi \eta kT}\frac{Q}{d}$$
(10)
Here, $`m`$ is the molecular weight of the gas, $`k`$ is Boltzmann’s constant, and $`T`$ is the absolute temperature. Low molecular weight is desirable, making helium a good choice for the transport gas. To get maximum laminar flow rate for a given capillary diameter $`d`$, Equation (10) is solved, using the estimated Reynold’s number of 1000:

$$Q_{max}(d)=\frac{1000d\pi kT\eta }{4m}$$
(11)
Taking $`Q=Q_{max}(d)`$ minimizes the transport time. Using the flow rate $`Q_{max}(d)`$ as a function of $`d`$, the head pressure $`P_p`$ needed to drive the gas flow can be determined from Equation (5). The efficiencies for each component of the system are computed for a given value of $`d`$ in order to find the optimum diameter, as shown in Figure 4. This figure illustrates that the optimal transport capillary diameter is approximately 3.25 mm if the decay chamber is at 125 kPa(absolute), giving an approximate total efficiency of $``$ 2.5%. Higher efficiencies may be achievable if the decay chamber pressure can be lowered using an exhaust pump, but this has not been tested in SNO.

The final system configuration, described in Table 1, provided a workable compromise between overall efficiency and robust operation.

## 4 Decay Chamber Design

The decay chamber is designed to fulfill three criteria: it must (a) be capable of detecting 8Li decays with a high efficiency for event identification during analysis, (b) provide the best possible electron spectrum to the SNO detector, and (c) it must be clean, safe, and robust to protect the ultrapure SNO detector media. These last two criteria are somewhat in competition, since the best electron spectrum requires a thin-walled chamber to minimize energy loss of the escaping beta, while safety demands that the wall be thick enough to resist buckling when the chamber is under pressure underwater (at a maximum depth of 18 meters of D2O). The main volume of the decay chamber is made of spun stainless steel in a sphere, which provides excellent resistance to compressive forces, good cleanliness for exposure to the water, and elastic properties such that a mechanical failure will tend to deform but not rupture the container. The sphere has a radius of 6.35 cm, with a wall thickness of 0.06 cm, corresponding to an electron energy loss of approximately $``$1 MeV for the electron energies of interest.

A diagram of the decay chamber is shown in Figure 5. Above the spherical decay volume is a conical ‘neck’ which contains the in- and out-flow gas ports that are connected to the umbilical (described in Section 3). Also in this neck is a thin acrylic window atop which sits a 1.5” Hamamatsu R580 photomultiplier tube (PMT) viewing the decay volume. The PMT, voltage supply, and gas fixtures are contained inside a sealed stainless steel endcap sitting atop the neck.

Event identification is performed by utilizing the scintillating properties of the helium transport gas. The 2$`\alpha `$ decay following the 8Li beta decay creates two alpha tracks nearly back-to-back with a mean energy of 1.6 MeV per particle. In the helium gas, each of these particles has a range of about 4 cm, so a large proportion of the alpha energy is deposited in the gas.

The scintillation signal is not large  and consists mostly of UV light which is outside the acceptance of the PMT photo-cathode. Moreover, the PMT photo-cathode coverage is a small fraction of the decay volume, so measures were taken to enhance the signal. The inside of the sphere is painted with a reflective white titanium oxide based paint (Bicron BC620), and has a thin coating of tetraphenyl butadiene (TPB) which acts as a wavelength shifter. In addition, a small admixture of N2 gas ($`0.1\%`$) is added to the helium gas stream to act as an additional wavelength shifter.

The scintillation light from helium is known to be prompt  and the TPB and N2 wavelength shifters are also fast, so the scintillation signal was fast enough to be included in the 400 ns coincidence window in the SNO electronics. The tagging PMT signal pulse height was digitized and recorded in a spare SNO electronics channel.

Decay tagging by gas scintillation has several advantages over other detection methods (e.g. proportional counters, ZnS scintillation). The decay volume need not assume a specific shape, which allows construction of a thin-walled design. (Cylindrical shapes, such as those required for proportional counters, require a larger wall thickness for the same volume container.) Gas scintillation is much faster than many other methods, allowing the decay tag to be measured in coincidence with the associated (SNO detector) event, and making subsequent analysis simple and robust.

## 5 Source Performance And Application

The 8Li source has been successfully tested and deployed in the SNO detector. Calculations using the described gas model and data from short-transfer line bench tests yield a prediction for the final 8Li yield of 0.75 $`\pm `$ 0.25 decays per second. A yield of approximately 0.5 decays per second was observed with the source deployed within the D2O volume of the SNO detector. An unexpected feature of the system behaviour was that the source took $``$2 hours to climb to the maximum yield. This is likely due to the salt temperature reaching equilibrium, increasing the pickup efficiency of the aerosol by the helium gas stream. This equilibration was not observed in short run bench tests. The yield provided by the source allowed a collection of some 50 000 tagged 8Li events to be observed.

Bench tests, using coincidences of external scintillators to identify 8Li electrons, allowed a measurement of the tagging efficiency of the decay chamber. This efficiency was estimated at 92$`\pm `$5% of all 8Li decays in the ball. This tagging efficiency remained quite robust with changes in decay chamber pressure. The tag signal was prompt ($``$150 ns after the decay), sufficient to record the signal within the 400 ns coincidence window of the SNO electronics. The PMT tagging signal from each 8Li decay was approximately 10-20 photoelectrons.

The 8Li system generated signals created by contaminant radioisotopes. Evidence of these contaminants was indicated in bench tests, and was confirmed during deployment in the SNO detector. In addition to 0.5 $`s^1`$ 8Li decays, approximately 8 $`s^1`$ of $``$5 MeV events were observed. Analysis of the signal indicated that the majority of this contamination was likely to be 16N, an isotope that decays by beta-gamma decay with a 6.13 MeV gamma and a half-life of 7.13 s. This isotope was likely produced by neutron reactions such as 19F (n,$`\alpha `$) and 16O (n,p) on Teflon and plastic components in the target chamber, and then carried to the decay chamber by the gas stream.

This contamination signal, however, did not interfere with the 8Li signal, because the majority of these events generated no tag pulse. The beta-gamma contaminants have only a small chance ($``$0.6%) for the beta particle to interact with the acrylic window at the top of the decay volume, generating 1-2 photoelectrons in the PMT by Cherenkov radiation. (The gas scintillation signal created by gammas or betas is much too small to be detected.) This contamination was reduced even further by using the time and integrated charge of the tagging signal to discriminate against the contaminant events. The PMT anode pulses are much smaller than those due to alpha particles from 8Li decays, allowing a clear separation of signals. In addition, the PMT timing is different (i.e. has a faster rise time) for $`\beta `$ detection in the chamber. These properties were tested by running CO2 gas through the system, generating large quantities of 16N isotope directly (in a manner completely identical with that described in Reference ). This test showed that the decay chamber detected 16N in a fashion consistent with the observed contaminant activity.

Event discrimination is illustrated in Figure 6, a scatter plot of PMT time and charge for the 8Li runs in the SNO detector. The $`\beta `$-tagged contaminant events are clearly distinguishable from the 8Li $`\alpha `$ signals. Figure 7 shows the Nhit distribution for position reconstructed events of all 8Li runs. It also shows Nhit distributions for the remaining events with a tag requirement and additional software cuts on PMT charges and times to distinguish 8Li from $`\beta `$-tagged contaminant events such as 16N. Using event discrimination, we estimate that the total probability of a single contaminant event to be misidentified as an 8Li event to be less than 0.005%. With the observed event rates, we conservatively estimate that events due to contaminants make up less than 0.1% of the final 8Li data set.

The Nhit response of the SNO detector to tagged 8Li decays is shown in Figure 8 for a run with the source situated at the center of the SNO detector. The Nhit distribution is shown to agree well with a simulation of the source using the SNO Monte Carlo code and an energy scale established with the 16N source . This comparison forms one of the suite of calibrations used to test energy response in the full SNO analyses . Important systematic uncertainties associated with such a comparison are uncertainties in the $`\beta `$ energy loss through the stainless steel wall and optical properties (shadowing and reflectivity) of the chamber. These uncertainties are on the order of 2% for the energy scale, with further investigation underway. An intriguing application of these data is direct comparison with the observed SNO 8B neutrino spectrum. Such a comparison has significant potential  to reduce systematic uncertainties associated with 8Be states common to both the 8B and 8Li decay channels.

The source has been deployed in multiple locations within the SNO D2O volume, and the resulting signals have been reconstructed in the same manner as neutrino data. A sample of reconstructed event positions is shown in Figure 9. Comparison of the source position and reconstructed event position has been used to evaluate the accuracy and resolution of position reconstruction algorithms .

Finally, the 8Li data comprise one of the tests used to evaluate data reduction techniques employed to remove instrumental background from the SNO detector data. The set of tagged 8Li events is used as a representative sample of data with distributions nearly identical to that of neutrino interactions in the detector. By counting the fraction of 8Li events removed by the data reduction cuts, the number of neutrino events falsely sacrificed by the cuts may be estimated, as described in Reference .

## 6 Conclusions

The 8Li source developed for SNO provides several necessary functions required to test the SNO detector response. Despite difficulties in producing and delivering the short-lived isotope, careful optimization allows an appreciable quantity of 8Li to be safely delivered over moderate distances.

Helium scintillation provides a mechanism to tag 8Li decays in a robust fashion. This method provides a signal that is sufficiently prompt to do fast coincidence with other signals, and has a high tagging efficiency. This signal can be used to discriminate against beta-gamma radioactive contaminants produced by the neutron source.

The source has been used successfully to evaluate SNO detector responses for energy, vertex reconstruction, and event topology . These calibrations are an important part of the SNO calibration programme .

## 7 Acknowledgements

The authors wish to thank R. Deal, P. Dmytrenko, J. Fox, E. Gaudette, V. Koslowsky, and E. Hagberg for their expertise and assistance in construction of the target and transfer systems, K. McFarlane, P. Liimatainen and D. Sinclair for their hard work in construction of the 8Li decay chamber, and C. Hearns and F. Duncan, for their work in deployment. We thank P. Huffman for kindly providing the TPB wavelength shifter. We are indebted to J.R. Leslie for suggesting the 11B(n,$`\alpha `$) reaction as a production method, and E. Bonvin for first pointing out the possiblilites of using 8Li as a calibration source. Many thanks go to the entire SNO collaboration for their support and help in making this project possible.

The SNO project has been supported in Canada by the Natural Sciences and Engineering Research Council, Industry Canada, National Research Council of Canada, Northern Ontario Heritage Fund Corporation and the Province of Ontario, in the United States by the Department of Energy, and in the United Kingdom by the Science and Engineering Research Council and the Particle Physics and Astronomy Research Council. The heavy water has been loaned by AECL with the cooperation of Ontario Power Generation. The provision by INCO of an underground site is greatly appreciated.

Table 1: Physical parameters of 8Li gas transport system. Ranges indicate approximate working values during operation, and are in rough agreement with the gas model given in the text.

Figure 1: The 8Li and 8B decay schemes. Data taken from Reference .

Figure 2: The basic features of the 8Li calibration source. Radioisotope is produced using fast neutrons impingent on a 11B-enriched target. Helium gas containing NaCl aerosol is used as a transport medium to carry 8Li from the target chamber to the decay chamber, where 8Li decays give off $`\beta `$-particles, which escape into the SNO detector volume, and $`\alpha `$-particles, which are tagged within the decay chamber.

Figure 3: The target (production) chamber, in situ about the neutron source. Gas flows into the volume through holes in perforated distribution plates. A thin cylinder of 11B-enriched boron in polyethylene sits inside the gas volume, and a thin layer of boron nitride coats the ouside of the gas volume. 8Li atoms recoil into the gas volume and are swept by laminar flow to a series of holes around the circumference of the chamber, which connect to the output line.

Figure 4: Partial efficiencies for exiting the target chamber, transport across the capillaries, and the efficiency for the decay in the trigger chamber are shown for different transport tube diameters. Also given is the overall efficiency.

Figure 5: The 8Li decay chamber. Two thin hemispheres of spun stainless steel are electron-welded together to create the decay volume at the bottom. These are in turn welded to the “neck” which holds the gas fixtures, and has a central hole for a 4.86 cm diameter acrylic window. Atop this window is a 3.81 cm diameter PMT. The PMT and gas fittings are contained within a sealed stainless steel endcap. Gas flows both in and out near the top of the volume. The top of this endcap holds the end of the umbilical (sealed with compressed O-rings). The source is supported from above by the SNO source manipulator system. Measurements are given in centimeters.

Figure 6: The integrated charge and time of signals coming from the decay chamber PMT for 47823 tagged events. Contaminant $`\beta `$-particles are clearly distinguishable from the $`\alpha `$ scintillation signal from 8Li. The line shows the cut used to discriminate against contamination events.

Figure 7: The Nhit spectra are shown for all source-related events within the detector during 8Li runs. (Note this spectrum is dominated by the 16N gamma contamination described in the text.) Also shown are spectra after cuts on the decay chamber PMT signal designed to distinguish tagged contaminant events from the tagged 8Li events.

Figure 8: The Nhit spectrum of tagged 8Li events is shown compared to a SNO Monte Carlo prediction of the source. The Monte Carlo energy scale was tuned using the 16N calibration source.

Figure 9: The left plot shows the x-z projection of reconstructed positions of $``$15k events for two placements of the 8Li source in the SNO detector, (x=-16 cm,y=24 cm,z=0 cm) and (x=-16 cm,y=24 cm,z=450 cm). The right plot shows a histogram of the distance between the source center and the reconstructed event positions for the z=0 cm position (solid circles) and the z=+450 cm position (open circles). (The center of the acrylic vessel is defined to be at (0,0,0) and the z-axis is vertically upwards.)


text:
# Sampling Chaotic Trajectories Quickly in Parallel

## I Introduction

Chaotic systems cannot be predicted for very long times because of the exponential divergence of nearby trajectories. Associated with the divergence of trajectories is a lack of history dependence; the current behavior of the system is not dependent on the past behavior. The absence of history dependence can be understood in various ways. Here I take a computational perspective on chaotic systems and analyze trajectories in terms of computational complexity. If, using a massively parallel computer, a typical long trajectory can be manufactured in far fewer parallel steps than the actual length of the trajectory then the trajectory lacks history dependence or historical complexity. Alternatively, if parallelism does not allow one to generate a typical trajectory much more quickly than its actual length then the trajectory displays a complex history dependence.

These considerations are illustrated using the one-dimensional quadratic map,

$$x_{n+1}=rx_n(1x_n)f(x_n)$$
(1)
with $`x_n[0,1]`$ and $`0<r<4`$. The objective is to produce typical trajectories for the map but since computing devices are necessarily restricted to finite precision we are really interested in generating typical pseudotrajectories NuYo . A pseudotrajectory of accuracy $`\delta `$ is a sequence $`\{y_n|n=0,\mathrm{},t\}`$ such that for all $`n`$ ($`0n<t`$),

$$|y_{n+1}f(y_n)|<\delta .$$
(2)
We suppose that the parallel computer works to a precision substantially better than $`\delta `$ so that bounds such as given in Eq. 2 can be checked with reasonable certainty.

For chaotic dynamics with a positive Lyapunov exponent $`\lambda `$, a pseudotrajectory and an exact trajectory that are initially equal remain close only for a time that is roughly given by $`(\mathrm{log}\delta )/\lambda `$. Nonetheless, for small $`\delta `$ a pseudotrajectory will have nearly the same statistical properties as a real trajectory and, in any case, numerical results about chaotic systems are learned from pseudotrajectories not real trajectories.

The goal then is to produce a pseudotrajectory chosen from the uniform distribution over pseudotrajectories. A typical pseudotrajectory of length $`t`$ can be generated using one processor in time linear in $`t`$ by iterating the map using arithmetic of precision much better than $`\delta `$ and then adding a noise term on each step chosen from the uniform distribution on $`[\delta ,\delta ]`$. Could a typical trajectory be produced in far fewer than $`t`$ parallel steps with the aid of many processor?

The model of parallel computation implicit in this question is the PRAM (parallel random access machine), the standard model in the theory of parallel computational complexity GrHoRu . A PRAM is an idealized, fully scalable device with many identical (except for distinct integer labels) processors. Processor all run the same program and all communicate with a global memory in unit time. Massive parallelism is envisioned here however the number of processors is required to be polynomially bounded in $`t\mathrm{log}(1/\delta )`$, the effective number of (binary) degrees of freedom of a pseudotrajectory of length $`t`$ and accuracy $`\delta `$.

In the next section we show how to produce a typical trajectory in parallel using Monte Carlo path sampling. The procedure is correct but inefficient. We then describe how simulated annealing together with path sampling can produce a typical pseudotrajectory in parallel time that scales linearly in $`\mathrm{log}t`$ and polynomially in $`\mathrm{log}(1/\delta )`$ using a number of processors that is polynomial in $`t\mathrm{log}(1/\delta )`$.

## II Path sampling of pseudotrajectories

The method described in this section is based on path sampling ideas put forward by Chandler and collaborators DeBoCsCh . The uniform probability density for pseudotrajectories , $`𝒫(y_0,y_1,\mathrm{},y_t)`$ is given by

$$𝒫(y_0,y_1,\mathrm{},y_t)=p(y_0)\underset{n=1}{\overset{t}{}}P(y_n|y_{n1})$$
(3)
where

$$P(y^{}|y)=\{\begin{array}{cc}\frac{1}{2\delta }\hfill & \text{if }|y^{}f(y)|<\delta \hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(4)
and $`p(y)`$ is the invariant distribution.

A simple Monte Carlo procedure can be used to sample paths from $`𝒫`$. Given an existing pseudotrajectory, a single time $`n>0`$ is chosen and a proposal for a new value for $`y_n`$ is obtained according to

$$y_n^{}=f(y_{n1})+ϵ$$
(5)
where $`ϵ`$ is chosen from the uniform distribution on $`[\delta ,\delta ]`$. The proposed value is accepted if it is also the case that $`|y_{n+1}f(y_n^{})|<\delta `$. It is straightfoward to verify that this Monte Carlo procedure satisfies detailed balance with respect to $`𝒫`$. The question of ergodicity of the Markov chain in the space pseudotrajectories is less clear. However, even if ergodicity holds, the actual mixing time for the Monte Carlo procedure would be long when $`\delta `$ is small since the time to obtain an independent pseudotrajectory is at least as great as $`1/\delta ^2`$, the time to diffuse a distance order one given a step size of $`\delta `$. To be considered an efficient process for generating pseudotrajectories , the parallel time should increase no more rapidly than some power of the logarithm of $`1/\delta `$. This goal can be achieved by first using a simulated annealing procedure to produce a pseudotrajectory that is random on long time scales then to use the above Monte Carlo path sampling to randomize the small scales.

## III Simulated annealing for pseudotrajectories

Long pseudotrajectories are constructed by independently generating many segments or short trajectories and then “welding” the segments together. In the welding step the discontinuity between successive segments is eliminated by simulated annealing. In addition to simulated annealing, it is sometimes necessary to extend segments to obtain a weld to the next segment.

The fundamental time scale in the system is

$$\tau =(\mathrm{log}\delta )/|\lambda |,$$
(6)
the time required for typical errors to grow to be order one. The length of the segments, $`K`$ used in the construction should be longer than the fundamental time $`\tau `$ so that the beginning and end of each segment is uncorrelated. We also need to allow for extensions at either end of the segment and for a “warm-up” so that the initial point of the segment is chosen from invariant distribution. Thus, in practice, for each segment we start with a random number and iterate the map $`L>K`$ times, choosing the segment of length $`K`$ from a predetermined part of the longer sequence of length $`L`$.

Having made a collection of segments, we now attempt to weld them together into a long pseudotrajectory . This is done in such a way that the initial value or anchor point of each segment is held fixed. The discontinuity between successive segments is annealed until all errors are less than $`\delta `$. The Monte Carlo annealing procedure is designed to lower the error $`e(y_{n1},y_n,y_{n+1})`$ associated with three successive elements, $`y_{n1},y_n`$ and $`y_{n+1}`$,

$$e(y_{n1},y_n,y_{n+1})=|y_nf(y_{n1})|+|y_{n+1}f(y_n)|.$$
(7)
If $`e(y_{n1},y_n,y_{n+1})<\delta `$ for every $`n`$ then we have a pseudotrajectory. For each time $`n`$, with the exception of the anchor points, the Monte Carlo annealing procedure begins with a measurement of $`e=e(y_{n1},y_n,y_{n+1})`$. If $`e<\delta `$ nothing is done. Otherwise a new value $`y_n^{}`$ is proposed,

$$y_n^{}=y_n+\stackrel{~}{ϵ}$$
(8)
where $`\stackrel{~}{ϵ}`$ is chosen as a Gaussian random variable with mean zero and standard deviation $`e/2`$. If $`e^{}=e(y_{n1},y_n^{},y_{n+1})`$, is less than $`e`$, the proposal is accepted as the new value for $`y_n`$. If the error increases, the proposal is accepted with probability $`e^{\beta (e^{}e)}`$. The value of the inverse temperature for each Monte Carlo step and is taken to be $`\beta =e/2`$ so that the acceptance ratio is independent of the size of the error. In a single Monte Carlo sweep, all except the initial and final values of each segment are processed using the above procedure. Given $`t/2`$ processors this can be done in constant parallel time by first processing the even and then the odd values of the time $`n`$.

In the chaotic regime ($`\lambda >0`$) the annealing procedure should yield a valid pseudotrajectory for large enough $`K`$ and sufficiently many Monte Carlo sweeps. In practice, however, some welds require a very large number of sweeps. Specifically, the probability distribution for the number of sweeps needed to achieve a weld has a long tail leading to parallel running times for creating pseudotrajectories that are dominated by the few most difficult welds. Two additional kinds of steps, called forward shifts and backward shifts cure this difficulty. Suppose segment $`s`$ together with its final condition, the first element of segment $`s+1`$, is not fully annealed after a predetermined number of annealing sweeps. Then segment $`s`$ is restored to its original state and either it is extended forward one step or segment $`s+1`$ is extended backwards one step. In the case of a backwards shift, the element prepended to segment $`s+1`$ is considered a new anchor points and serves as the new final condition for segment $`s`$. The net effect of either a forward or backward shift is that the discontinuity between segments $`s`$ and $`s+1`$ occurs with a different pair of numbers. Annealing sweeps and shifts are interleaved, a fixed number of Monte Carlo annealing sweeps are attempted and if all errors are not less than $`\delta `$, a shift is done. The process is repeated until a satisfactory weld is achieved. Successive shifts are alternately of the forward and backward type. In Sec. V, I show that the combination of Monte Carlo annealing and shifts produces a pseudotrajectory in $`O(\mathrm{log}t)`$ parallel steps.

Shifts serve several purpose. First, they simply provide for the possibility of more Monte Carlo sweeps, though if this were their only function it could be accomplished by directly increasing the number of sweeps. Second, shifts permit the algorithm to perform properly for periodic orbits or nearly periodic stretches of aperiodic orbits. The annealing procedure by itself cannot generate long pseudotrajectories for periodic orbits since welds are often attempted between segments that are out of phase with one another. Adding shifts to the annealing procedure insures that periodic pseudotrajectories will be correctly generated. For example, consider the case of a period two orbit, a single forward or backward shift of some segments will insure that all welds are satisfactory. For period $`d`$ orbits, as many as $`d1`$ shifts are necessary to insure that all welds are satisfactory.

Shifts may also provide padding around hard to weld regions of a trajectory. During a shift, new points are added to one end of a segment but no points are removed. Thus shifts do not bias the pseudotrajectory against difficult to weld regions in the invariant measure. For example, it is observed that if the final condition for a segment is very near the maximum of the support in the invariant measure at $`r/4`$ then one or more backward shifts are usually necessary so that the point near $`r/4`$ is surrounded by a region of small errors and is not involved in the annealing process.

## IV Full parallel algorithm for pseudotrajectories

This section provides the details of the parallel algorithm for producing typical pseudotrajectories that combines the path sampling method of Sec. II and the annealing procedure of Sec. III. First the annealing procedure is used to generate a pseudotrajectory and then path sampling is used to further randomize it. The algorithm is controlled by several parameters: $`t`$ is the total length of the desired pseudotrajectory , $`\delta `$ is the desired accuracy, $`K`$ is the length of each segment, $`K^{}`$ is the warm-up length, $`E`$, an even number, is the maximum number of shifts that are attempted, $`M_1`$ is the number of Monte Carlo annealing sweeps carried out between shifts and $`M_2`$ is the number of Monte Carlo path sampling sweeps. The algorithm is described below:

1. In parallel, generate $`S=t/K`$ sequences $`\{x_m^{(s)}\}`$ each of length $`L=K^{}+K+E`$, with $`s=0,\mathrm{},S1`$ and $`m=0,\mathrm{},L1`$. The initial value of each sequence is a uniform random number on $`(0,1)`$ and subsequent values are obtained by iterating the map $`L1`$ times to precision much greater than $`\delta `$. This step requires $`O(L)`$ parallel time.
2. These $`S`$ sequences are used to define $`ES`$ segments $`\{y_n^{(s,q)}\}`$ each of length $`K`$, where the index $`q`$, $`0q<E`$ gives the number of shifts. $`y_n^{(s,q)}=x_{n+E/2+K^{}+q/2}^s`$ for $`0n<K1`$ while $`y_{K1}^{(s,q)}=x_{E/2+K^{}q/2+K^{}}^{s+1}`$. Note that the final point in segment $`y^{(s,q)}`$ is a taken from sequence $`s+1`$.
3. In parallel, for each $`s<S`$, and each $`q<E`$, anneal segment $`y^{(s,q)}`$. The annealing procedure consists of $`M_1`$ Monte Carlo sweeps. During a single annealing sweep first the even and then the odd elements of the segment are updated in parallel. The initial and final points $`y_0^{(s,q)}`$ and $`y_{K1}^{(s,q)}`$ are held fixed during the annealing procedure. A single Monte Carlo update of the point $`y_m^{(s,q)}`$ consists of the following procedure:

* Compute $`e=e(y_{m1}^{(s,q)},y_m^{(s,q)},y_{m+1}^{(s,q)})`$ from Eq. 7.
    * If $`e<\delta `$, do nothing. If $`e>\delta `$ propose a new value $`y^{}=y_m^{(s,q)}+\stackrel{~}{ϵ}`$ where $`\stackrel{~}{ϵ}`$ is a Gaussian random variable with mean zero and standard deviation $`e/2`$.
    * Compute $`e^{}=e(y_{m1}^{(s,q)},y^{},y_{m+1}^{(s,q)})`$. If $`e^{}e`$ accept the proposed change, $`y_m^{(s,q)}y^{}`$. If $`e^{}>e`$ accept the proposed move with probability $`\mathrm{exp}[\beta (e^{}e)]`$ where $`\beta =e/2`$.

The parallel time required for this step is $`O(M_1)`$.
4. In parallel, for each $`s<S`$ find $`Q(s)`$, that smallest value of $`q`$ such that the segment is successfully annealed. The annealing is successful for this segment if, for all $`0<m<K1`$, the errors are sufficiently small, $`e(y_{m1}^{(s,q)},y_m^{(s,q)},y_{m+1}^{(s,q)})<\delta `$. If for any $`s`$, annealing is unsuccessful for all $`qE`$, the algorithm fails. This step can be carried out in constant parallel time.
5. The full pseudotrajectory $`y^{}`$ is a concatenation of sequences obtained from each original sequence $`x^{(s)}`$. The contribution to the pseudotrajectory from $`x^{(s)}`$ is the concatenation of $`\{x_m^{(s)}|E/2+K^{}Q(s1)/2mE/2+K^{}+Q(s)/2\}`$ and $`\{y_m^{(s,Q(s))}|0<m<K1\}`$. The first of these sequences is composed of the anchor points and the second sequence is the annealed segment. To obtain a pseudotrajectory of length exactly $`t`$, the pseudotrajectory obtained above is simply truncated after $`t`$ steps.
6. The path sampling Monte Carlo procedure described in Sec. II further randomizes $`y_m^{}`$. During a single sweep, first all the even and then all the odd elements of the pseudotrajectory are updated in parallel. The number of sweeps is $`M_2`$. On each Monte Carlo step a new value for $`y_m^{}`$ is proposed according to Eq. 5 and accepted only if the trajectory is still a pseudotrajectory within error $`\delta `$. The randomization step requires parallel time $`O(M_2)`$.

## V Validity and complexity of the parallel algorithm

The central questions addressed in this section are (1) whether the algorithm succeeds in creating a pseudotrajectory, (2) how the scaling of the number of parallel steps depends on the length and accuracy of the pseudotrajectory and the parameter $`r`$ of the map and (3) whether the algorithm samples the uniform distribution on pseudotrajectories. A sequential algorithm that carries out the annealing and path sampling routines one segment at a time was used to study these questions. In the simulations reported below, the parameters are chosen to be $`M_1=M_2=5\tau ^2`$, $`K=5\tau `$ and $`K^{}=1000`$. The assumption behind these choices is that memory is lost on a time scale $`\tau `$ so that placing independently chosen anchor points separated by $`K=5\tau `$ is satisfactory. The annealing process that welds successive segments is expected to influence a region whose length is order $`\tau `$. Since information is transmitted diffusively by local Monte Carlo moves, having the number of Monte Carlo sweeps scale as $`\tau ^2`$ should suffice.

The annealing stage of the parallel algorithm can be studied one segment at a time since each segment is independently annealed. First, I observed that, given enough shifts, the annealing step always produced a successful weld. The choice of the maximum number of shifts $`E`$ for the annealing stage must be large enough to make the failure probability for the whole algorithm small. For long trajectories, $`t\tau `$, the choice of $`E`$ is determined by the tail of the distribution of the number of shifts, $`Q`$ required to obtain a weld since the whole procedure fails if even one segment is not successfully annealed. Suppose $`C()`$ is the cumulative probability distribution for $`Q`$. An estimate of the maximum number of shifts, $`E`$ needed to insure all segments are annealed is given by the relation

$$(t/K)(1C(E))<1.$$
(9)
Figure 1 shows $`\mathrm{log}_{10}(1C(Q))`$ vs. $`Q`$. These data were collected for the case of $`r=3.7`$ (the period doubling transition to chaos occurs at $`r=3.5699\mathrm{}`$) and the six curves from left to right are for $`\delta =10^5`$ through $`10^{10}`$, respectively. Each curve is obtained from annealing $`10^5`$ segments except for the $`\delta =10^{10}`$ curve which is obtained from $`6\times 10^4`$ segments. For $`r=3.7`$ the Lyapunov exponent is $`\lambda =0.354`$ and so, for example, with $`\delta =10^7`$, $`\tau =45.5`$, $`K=228`$ and $`M_1=10366`$. Over a reasonable range following an initial transient and before the noise becomes large, the data falls on straight lines suggesting that the distribution is asymptotically exponential, $`C(Q)1\mathrm{exp}(Q/\sigma )`$. Equation 9 then implies that $`E\sigma \mathrm{ln}t`$ and we can conclude that the parallel running time is $`O(\mathrm{log}t)`$ since no other contribution to the running time depends on the overall length of the pseudotrajectory.

How does the decay constant $`\sigma `$ and thus the running time depend on the choice of the accuracy $`\delta `$. Figure 2 shows $`\sigma `$ vs. $`\delta `$ on a logarithmic scale for the case $`r=3.7`$ and suggests that $`\sigma `$ is a polynomial function of $`\mathrm{log}\delta `$. The other simulation parameters, $`L`$, $`M_1`$ and $`M_2`$ are also polynomial in $`\mathrm{log}\delta `$ so we conclude that the full algorithm has a running time that is polynomial in $`\mathrm{log}\delta `$ and linear in $`\mathrm{log}t`$.

I also considered two other parameter values for the quadratic map, $`r=3.6`$ where $`\lambda =0.183`$ and $`r=3.95`$ where $`\lambda =0.577`$. In both cases, the accuracy was set to $`\delta =10^7`$. The decay of $`(1C(Q))`$ appears to be exponential in both cases but with rather different values of the decay constant: $`\sigma =450,15`$ and $`9.5`$ for $`r=3.6,3.7`$ and $`3.95`$, respectively. Either $`\sigma `$ depends strongly on $`\lambda `$ or perhaps there are additional $`r`$ dependent factors controlling $`\sigma `$. For example, for $`r=3.7`$ and 3.8 the invariant measure has support on a single interval but for $`r=3.6`$ the support consists of two intervals.

The annealing stage of the algorithm creates a pseudotrajectory but it is not typical in the sense of being chosen from the distribution of Eq. 3. On long time scales, the pseudotrajectory is randomized by the random choice of initial conditions for each sequence. However, the individual errors, $`ϵ_n=y_{n+1}^{}f(y_n^{})`$ are not guaranteed to be independent random variables on the interval $`[\delta ,\delta ]`$. For example, anchor regions of the pseudotrajectory have errors much less than $`\delta `$. The hypothesis is that $`M_2=O(\tau ^2)`$ path sampling Monte Carlo sweeps are sufficient to randomize the short time scales and produce a typical pseudotrajectory from the pseudotrajectory produced by the annealing stage. To check this hypothesis, I computed mean values and autocorrelation functions for errors and cross correlations between errors and values of the pseudotrajectory . The quantities $`ϵ_n`$, ($`ϵ_n^2\delta ^2/3`$), $`ϵ_{n+1}ϵ_n`$, $`ϵ_{n+1}y_n^{}`$ and $`(ϵ_{n+1}^2ϵ_n^2\delta ^4/9)`$ were all found to be zero within error bars for the case $`r=3.7`$ and $`\delta =10^7`$. Here the angled brackets indicate an average over segments and over $`n`$. The vanishing of these quantities is a necessary but not sufficient condition that pseudotrajectory is chosen from the uniform distribution described by Eqs. 3 and 4. More work is needed to firmly establish that the algorithm with $`O(\tau ^2)`$ path sampling sweeps samples the uniform distribution to good approximation. However, even without the path sampling stage, the pseudotrajectories produced by the annealing stage are typical in a different sense. As shown in NuYo , any pseudotrajectory shadows an exact trajectory (i. e. remains close to over its entire length) though possibly for a larger value of $`r`$.

## VI Conclusions

I have exhibited a parallel algorithm that generates pseudotrajectories of the quadratic map. Numerical evidence suggests that the parallel time required to generate a typical pseudotrajectory increases linearly in $`\mathrm{log}t`$ and polynomially in $`\mathrm{log}(1/\delta )`$ though more work would be required to establish these scalings with certainty. Essentially the same parallel algorithm can be applied to other one-dimensional and higher dimensional maps. It would be interesting to explore whether the annealing/shift procedure is sufficient to efficiently sample pseudotrajectories for other maps.

Since there is little demand for very long pseudotrajectories of the quadratic map, the fast parallel algorithm is probably not of practical value. The significance of its existence and complexity is that it characterizes the history dependence of the map. The existence of a fast parallel simulation is a strong statement against history dependence since it shows that the logical path from independent random numbers (used to drive the Monte Carlo procedures) to a typical pseudotrajectory is much shorter than the length of the pseudotrajectory. The length of this logical path is one measure of the potential for generating historical complexity. In very few logical steps, very little complexity can be arise. The idea that complexity tends to emerge at the “edge of chaos” Lang  is born out here since the basic time scale $`\tau `$ for the parallel algorithm diverges when the Lyapunov exponent vanishes. An appealing feature of characterizing systems by computational complexity is that very different systems in statistical physics systems, for example diffusion limited aggregation MoMaGr97 , sandpiles MoNi  or the Bak-Sneppen model MaLi01 , can be compared to one another within the same framework.

###### Acknowledgements.

This work was supported in part by NSF Grant DMR 9978233.


text:
# Extended Comment on Language Trees and Zipping

## Appendix A Experimental Results

In this appendix, I briefly describe the experimental results comparing Naive Bayes to zipping. The data set used was the 18828 version, available from http://www.ai.mit.edu/~jrennie/20Newsgroups. This version has had duplicates and most headers removed. No additional processing of the data was done. For Naive Bayes, words were simply segmented at white-space boundaries using the perl “split” function.

For Naive Bayes, I computed for each topic the total number of times each “word” occurred in the topic. I used standard plus-1 smoothing, by, for each topic, adding one count for every word that occurred in any training set. The probability assigned to a given word was thus:

$$P(\text{word}|\text{topic})=\frac{\text{count}(\text{word in topic})+\{\begin{array}{cc}1\hfill & \text{if }\text{word}\text{ occurs in any topic}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}}{\text{total length of topic}+\text{total number of words that occur in any topic}}$$
(1)
The count of a word in a topic counts every occurrence of that word in the topic, counting each word multiple times if it appears multiple times in the topic. The total length of a topic counts every occurrence of every word in the topic, counting multiple times for words that occur multiple times. The total number of words that occur in any topic counts each word that ever appears exactly once, no matter how many times it appears.

The probability of a document of length $`n`$ given the topic is simply

$$P(\text{document}|\text{topic})=\underset{i=1}{\overset{n}{}}P(\text{word}_i|\text{topic})$$

If a word occurs more than once in the document, its probability is multiplied in more than once. (This estimate neglects the prior on document length.)

We wish to find the most probable topic, given the document. The probability of a topic given a document is found through a simple application of Bayes’ law:

$`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}P(\text{topic}|\text{document})`$ $`=`$ $`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}{\displaystyle \frac{P(\text{document}|\text{topic})\times P(\text{topic})}{P(\text{document})}}`$
$`=`$ $`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}P(\text{document}|\text{topic})\times P(\text{topic})`$
where we can remove the term $`P(\text{document})`$ because it is independent of the topic. Furthermore, we assume a uniform prior on topics (which is fine for 20 newsgroups, but may be inappropriate for other tasks.) Thus, we get

$`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}P(\text{topic}|\text{document})`$ $`=`$ $`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}P(\text{document}|\text{topic})`$
$`=`$ $`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}\mathrm{log}P(\text{document}|\text{topic})`$
$`=`$ $`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}\mathrm{log}{\displaystyle \underset{i=1}{\overset{n}{}}}P(\text{word}_i|\text{topic})`$
$`=`$ $`\mathrm{arg}\underset{\text{topic}}{\mathrm{max}}{\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{log}P(\text{word}_i|\text{topic})`$
Thus, in practice, the algorithm is extremely simple. For a given document, for each topic, we simply apply Equation 1 to each word in the document; we take the log to avoid underflow; we compute the sum of those values. That gives us a value based on the log of the probability of the topic given the document. We find the topic that has the largest negative log probability (which corresponds to the highest probability.) Note that all of these techniques are well known and standard.

For zipping, I applied the algorithm of \[Benedetto et al.,2002\].

The code to split the data into test and training was approximately 70 lines. The code to implement Naive Bayes was also approximately 70 lines. No attempt to keep the code small was made. The code to implement the zipping algorithm was about 50 lines. Note that all of these file sizes are trivial. For instance, the code I usually use to do experiments of this nature is over 7000 lines (all written by me.) It took very roughly an hour to write, debug, and verify the correctness of each program. (The zipping program was actually harder to write, despite its smaller size, because it ran so slowly that it was more time consuming to debug and verify correctness, and because some additional effort was made to speed it up, resulting in roughly a factor of two increase in speed.)

I used the same training and test set for both experiments: I used 99% of the data for training, and 1% for testing. This left 192 test documents, which is easily large enough to reliably detect errors of the magnitude seen here. (I used a relatively small test set because the zipping algorithm was so slow.) The Naive Bayes algorithm made 26 mistakes on this test set, a 14% error rate, while zipping made 91 errors, a 47% error rate. The Naive Bayes algorithm required 5 minutes to run, elapsed time, while zipping required one hour and 25 minutes, elapsed time, on the same machine with no other programs in use. This was despite minor attempts to optimize the zipping algorithm and no attempts to optimize the Naive Bayes algorithm. Note that the Naive Bayes algorithm was implemented in perl, which is typically much slower than, say, C or C++. The zipping script was also written in perl, but most of the time was spent in gzip, which is written in C.

Thus, it appears that with approximately equal effort, and using well known techniques, one can perform text categorization over 3 times more accurately and 17 times faster by using Naive Bayes rather than the zipping technique.