text:
# 1 Introduction

## 1 Introduction

The large majority of the precision detectors of the muon spectrometer of the ATLAS experiment at the Large Hadron Collider (LHC) will be Monitored Drift Tube (MDT) chambers. These chambers consist of 2 multilayers, each built of 3 or 4 layers of densely packed drift tubes mounted on a support frame made of aluminum (compare figure 1). The drift tubes — which are also made of aluminum — have an outer diameter of 3 cm and a wall thickness of 400 $`\mu `$m; in the middle of each tube a gold-plated W-Re anode wire of 50 $`\mu `$m diameter is stretched. The anode wire is positioned only at the two tube ends by precision end-plugs; the sag of each chamber can be adjusted to follow the sag of the wire along the tube. The drift tubes are operated at a pressure of 3 bar absolute with an Ar:CO2=93:7 gas mixture at a gas gain of $`2\times 10^4`$. The area covered by a complete chamber ranges from $`1\mathrm{m}^2`$ to $`11\mathrm{m}^2`$. A total of 1194 MDT chambers will be used in ATLAS .

In collaboration with the Max Planck Institute (MPI) for Physics, Munich University (LMU) is responsible for the construction of 88 MDT chambers of type BOS (Barrel Outer Small — the second largest type used in the ATLAS barrel region) which consists of 432 drift tubes, arranged in $`2\times 3`$ layers of 72 tubes, and has a size111Some chambers will be narrower or have cut-outs — realized with shorter drift tubes — to accommodate support structures of the ATLAS detector or feedthroughs to inner components. of $`3.9\mathrm{m}\times 2.2\mathrm{m}\times 0.5\mathrm{m}`$. The chambers are built with high mechanical precision at MPI  and are then commissioned and calibrated at the Cosmic Ray Measurement Facility of Munich University.

The main objectives of the Cosmic Ray Facility are to test the chambers, to map the positions of the tube layers and of individual anode wires within the chambers, and to gain experience in the operation of the chambers and the calibration procedures.

An important design goal of the ATLAS muon spectrometer is the standalone measurement capability of muon momenta up to 1 TeV with a relative error of less than 10%. This leads to the requirement that the anode wire positions within the MDT chambers must be known to within $`20\mu \mathrm{m}`$ .

## 2 The Cosmic Ray Measurement Facility Setup

A schematic view of the Cosmic Ray Measurement Facility is shown in figure 2. The whole setup is located in a climatized hall to ensure stable environmental conditions during the measurements. We define a right handed coordinate system, the $`x`$ coordinate in the direction along the wires, $`y`$ perpendicular to the wires in the chamber plane, and $`z`$ perpendicular to the chamber plane pointing upwards.

Three MDT chambers are operated simultaneously in the facility. Two of them, the so-called reference chambers, encompass the third chamber — which is to be tested — and provide track information for cosmic ray muons. The anode wire positions of the reference chambers have been mapped with a precision of $`2\mu \mathrm{m}`$ with an X-ray tomograph at CERN .

A coincidence of the three scintillator layers in the hodoscopes above and below the MDT chambers provides the trigger. The upper hodoscope consists of one layer of scintillator counters, the lower of two staggered layers displaced by half the width of a counter. Only muons with a momentum greater than 600 MeV pass the 34 cm iron absorber above the lower hodoscope and can generate a trigger signal. The trigger logic is divided into 5 segments along the $`x`$ coordinate to limit the inclination of the muons in the $`x`$-$`z`$ plane. The arrangement of the scintillator bars allows the measurement of muon tracks in the $`x`$ direction with a precision of 8 cm. The time of the muon transition is determined by the lower hodoscope with a precision of better than 800 ps. The average trigger rate is about 70 Hz.

The streamer tubes at the bottom of the setup are used together with the lower reference chamber to determine the multiple scattering angle in the $`y`$-$`z`$ plane in the iron absorber. This angle is used in combination with the angle between the track segments in the two reference chambers to estimate the momentum of the muon.

Two contact free alignment systems are used to continuously monitor the positions of all MDT chambers: The position of the reference chambers relative to each other is measured with an optical alignment system consisting of eight RasNik sensors , the position of the test chamber with respect to the upper MDT chamber with eight capacitive sensors . Both systems have a precision of better than $`5\mu \mathrm{m}`$ .

In addition to the chamber-to-chamber alignment, the internal chamber geometry is also monitored with integrated RasNik sensors (compare figure 1).

## 3 Drift Time Spectra

The position measurement in the MDT chambers is based on the measurement of the drift time of the electrons (from ionization by the incident muon along its trajectory) to the anode wire. The drift time is determined with respect to the external trigger system and translated to the impact radius (distance from the anode wire) via the so-called $`r`$-$`t`$ relation (see section 4); the drift tubes are only sensitive to the track coordinates in a plane perpendicular to the anode wire, not along it. The drift of the electrons depends on the gas composition, the gas density, and the electric field applied (i.e. the operating voltage), and is therefore characteristic of the gas mixture and the operating conditions.

After correcting for the signal propagation times along the anode wire and in the electronics and cables, the distributions of drift times are used as a first check that each tube is functioning as expected and that the detector response is homogeneous. For the simple case of uniform illumination of a tube, the density of hits at a given drift time is proportional to the drift velocity at the corresponding radius.

Events with more than one charged particle passing through the setup would lead to ambiguities in the hodoscope time measurement. Therefore, exactly one hit in the upper hodoscope and one pair of hits in overlapping counters in the lower hodoscope are demanded, which rejects 15 % of all triggers (see table 1).

The edges of the drift time spectra are parameterized  with the functions

$`F\left(t\right)=p_0+{\displaystyle \frac{A_0}{1+e^{\frac{t_0t}{T_0}}}}`$ $`\mathrm{for}\mathrm{the}\mathrm{rising}\mathrm{edge}\mathrm{and}`$ (1)
$`G\left(t\right)=p_m+{\displaystyle \frac{\alpha _mt+A_m}{1+e^{\frac{tt_m}{T_m}}}}`$ $`\mathrm{for}\mathrm{the}\mathrm{trailing}\mathrm{edge}\mathrm{of}\mathrm{the}\mathrm{spectrum},`$ (2)
as shown in figure 3. The parameters $`p_0`$ and $`p_m`$ correspond to the rates of accidental hits, $`A_0`$ to the height of the spectrum, and $`t_0`$ to the mid-time and $`T_0`$ to the steepness of the rising edge of the spectrum, which is determined by the resolution of the drift tube near the anode wire. Correspondingly, $`t_m`$ is the middle of the trailing edge of the spectrum and $`T_m`$ its steepness, which is determined by the decrease in pulse height near the tube wall222The small pulse height for muons passing a tube close to the tube wall is due to their short path inside the tube, which leads to smaller primary ionization.. The parameterization with $`\alpha _mt+A_m`$ accounts for the slope of the spectrum before its trailing edge that is due to the decrease of the drift velocity with increasing radius.

The maximum drift time, defined as $`\tau =(t_m+T_m)t_0`$, is characteristic of the drift properties of a tube. The distribution of the maximum drift time is shown in figure 4a for the tubes of one MDT chamber.

As part of the chamber quality control programme, it is checked that all tubes have the same drift time spectrum characteristics within errors, compare figure 4b.

## 4 The Relation between Drift Time and Drift Radius

The relation between the measured drift time and the corresponding drift radius ($`r`$-$`t`$ relation) is determined without the help of external spacial detectors and without assumptions on uniform tube illumination. After the check described in section 3, it is assumed that all tubes of an MDT chamber have the same $`r`$-$`t`$ relation.

Events with two hits at small drift times ($`|tt_0|<15\mathrm{ns}`$) in one MDT chamber are selected. By assuming that the muon went through the two corresponding anode wires and that the trajectory of the muon is a straight line, the impact radii in the other tubes of the chamber can be calculated. This procedure yields drift times for a discrete set of drift radii. The number of such points — typically 16 — is determined by the geometry of the chamber and the angular acceptance of the trigger. In addition to these points, the drift times corresponding to $`r=0`$ and to the inner tube radius are determined by the fit to the drift time spectra described in section 3. The $`r`$-$`t`$ relation is then approximated by a linear interpolation between the measured values.

## 5 Track Reconstruction

Within each MDT chamber a muon track segment is reconstructed from the measured hits. A track segment is assumed to be a straight line in the $`yz`$ plane333The third coordinate ($`x`$) is only used to correct for signal propagation times and for the wire sag and is taken from the hodoscope information.:

$$y(z)=m_nz+b_n,$$
(3)
where the index $`n`$ specifies the chamber. The parameters $`m_n`$ and $`b_n`$ are determined for each chamber by minimizing the $`\chi ^2`$ function

$$\chi ^2=\underset{i}{}\frac{\left(r_{track,i}r_{drift,i}\right)^2}{\sigma _i(r_{drift,i})^2},$$
(4)
where the index $`i`$ runs over all tubes contributing to the track segment, $`r_{track,i}`$ is the distance of the track segment to the wire, $`r_{drift,i}`$ is the measured drift radius, and $`\sigma _i`$ is the resolution of the tube for the given radius. At a test beam measurement using an external reference tracking system of silicon strip detectors the intrinsic single tube resolution was determined to be around $`100\mu \mathrm{m}`$ .

The tubes contributing to a track segment are selected by choosing one tube with a hit in each layer of the chamber, which is a sample of six tubes. From the impact radii in the two outermost tubes of this sample four possible track segment candidates are derived. If the difference between the distance of such a candidate to the wire and the measured drift radius does not exceed 2 mm in any of the other tubes, a track segment is fitted using this sample. If more than one track segment is found, the one with the smallest $`\chi ^2`$ is used. If no track segment is found, the procedure is repeated with subsamples of five tubes.

## 6 Chamber Alignment

The MDT chambers in the Cosmic Ray Facility are aligned mechanically with a precision of about $`100\mu `$m. A much more precise estimate of their relative positions is then obtained by comparing the muon track segments reconstructed in each chamber. The optical and capacitive alignment systems are used to verify that the setup does not move significantly during the data taking period. In the future, the alignment information will be included in the analysis to correct for time-dependent position shifts of the chambers.

The positions of the two reference chambers are measured relative to the test chamber.

* A shift $`v_y`$ of a reference chamber with respect to the test chamber results in a systematic shift between the track segments in the reference and test chamber and can be obtained from the average over all events,

$$v_y=b_{\mathrm{ref}}b_{\mathrm{test}},$$
(5)
where $`b_{\mathrm{ref}}`$ and $`b_{\mathrm{test}}`$ are the track segment parameters defined in equation (3). They are determined with a precision of 1 $`\mu `$m.
* A shift $`v_z`$ results in a relative shift of the track segments depending on the slope $`m`$. It is determined from a linear fit to the distribution of $`b_{\mathrm{ref}}b_{\mathrm{test}}`$ versus $`m`$ with an error of 10 $`\mu `$m.
* A shift $`v_x`$ along the $`x`$ coordinate can be ignored, as the MDT chambers are not sensitive to such a shift along the direction of their anode wires.
* The tilt angle $`\alpha `$ of a reference chamber around the $`x`$ axis is given by the systematic deviation of the two track segment slopes as

$$\alpha =m_{\mathrm{ref}}m_{\mathrm{test}}.$$
(6)
A precision of $`10^6`$ is achieved.
* To determine the angles $`\beta `$ and $`\gamma `$ corresponding to rotations of the reference chamber around the $`y`$ and $`z`$ axes, the $`x`$ coordinate measurement of the hodoscope is used. The shifts $`v_y`$ and $`v_z`$ are calculated for three one meter wide sections in $`x`$. From the dependences $`v_z(x)`$ and $`v_y(x)`$, the angles $`\beta `$ and $`\gamma `$ are obtained with precisions of $`8\times 10^6`$ and $`10^6`$, respectively.

In the further analysis, the parameters of track segments in the reference chambers are corrected according to the measured reference chamber positions.

## 7 Energy Estimation

The energy of cosmic muons traversing the Cosmic Ray Measurement Facility is estimated from multiple scattering — which depends on the muon momentum — in the chambers and the iron absorber.

If $`m_{\mathrm{ref},\mathrm{u}}`$ designates the slope of the track segment reconstructed in the upper reference chamber and $`m_{\mathrm{ref},\mathrm{l}}`$ the slope of the track segment reconstructed by the lower reference chamber, the distribution of $`\mathrm{\Delta }m=m_{\mathrm{ref},\mathrm{u}}m_{\mathrm{ref},\mathrm{l}}`$ will be the wider the lower the muon energy. The standard deviation $`\sigma _{\mathrm{\Delta }m}(E_\mu )`$ has been determined as a function of the muon energy by means of a Monte Carlo simulation and can be parameterized as

$`\sigma _{\mathrm{\Delta }m}(E_\mu )=\sigma _{\mathrm{}}+\sigma _0\left({\displaystyle \frac{600\mathrm{MeV}}{E_\mu }}\right)^\alpha ,`$ (7)
with $`\sigma _{\mathrm{}}=(2.0\pm 0.2)10^4`$, $`\sigma _0=(1.001\pm 0.008)10^2`$, and $`\alpha =1.04\pm 0.01`$.

Similarly, the deviation $`\mathrm{\Delta }y`$ of the position measured with the streamer tubes from the prediction obtained from the track segment in the lower reference chamber depends on the muon energy: For its width $`\sigma _{\mathrm{\Delta }y}(E_\mu )`$, a Monte Carlo simulation yields the values $`\sigma _{\mathrm{}}=(4.4\pm 0.1)`$ mm, $`\sigma _0=(183\pm 2)`$ mm, and $`\alpha =1.44\pm 0.01`$.

The probability density $`f`$ for measuring deviations $`\mathrm{\Delta }m`$ and $`\mathrm{\Delta }y`$ is then given by

$`f(E_\mu )={\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{1}{\sigma _{\mathrm{\Delta }m}(E_\mu )}}{\displaystyle \frac{1}{\sigma _{\mathrm{\Delta }y}(E_\mu )}}\mathrm{exp}\left({\displaystyle \frac{1}{2}}\left[\left({\displaystyle \frac{\mathrm{\Delta }m}{\sigma _{\mathrm{\Delta }m}(E_\mu )}}\right)^2+\left({\displaystyle \frac{\mathrm{\Delta }y}{\sigma _{\mathrm{\Delta }y}(E_\mu )}}\right)^2\right]\right).`$ (8)
For a given track, the value of $`E_\mu `$ that maximizes $`f(E_\mu )`$ is an estimate of the muon energy. This estimator is biased and has a limited resolution but is still useful for a selection of high momentum muons.

## 8 Measurement of Wire Positions

In this section, the method for the measurement of anode wire positions is presented. To determine wire positions, the muon track segments in the reference chambers are extrapolated into the test chamber, and the drift radii measured in the test chamber are compared with the track predictions.

The Cosmic Ray Measurement Facility is particularly sensitive to displacements $`\delta y`$ of the wires in the chamber plane. From this, a shift of entire tube layers in the $`y`$ direction and the mean spacing $`g`$ between the wires within a layer can also be determined.

Wire displacements $`\delta z`$ in the plane perpendicular to the chamber are accessible via tracks with different inclinations $`m`$ (cf. equation (3)). This measurement is less precise than the $`\delta y`$ determination owing to the limited angular acceptance of the trigger and the angular distribution of cosmic muons. During chamber production, $`z`$ displacement and tilts of entire tube layers are more difficult to control than individual wire displacements $`\delta z`$ with respect to the layer. Therefore, the $`z`$ displacement for each tube layer is also measured.

The measurements are performed at both ends of the chamber by selecting muons which passed the tested chamber within one meter of the corresponding end of the chamber using the hodoscope information. Because the wires are only supported at the tube ends and since the sag of the anode wires is known from wire tension measurements these two measurements completely determine the wire position.

### 8.1 Event Selection

Events with one reconstructed charged particle are selected as described in section 3. For the further analysis, the event selection proceeds in two steps: First, the reconstructed track segments are subjected to a set of cuts, and second, the individual hits in the tubes must pass certain quality criteria.

The presence of a reconstructed track segment with at least 5 hits is required in each of the three MDT chambers. The track segments must also roughly match, i.e., $`|b_{\mathrm{ref},\mathrm{u}}b_{\mathrm{ref},\mathrm{l}}|<4\mathrm{mm}`$ and $`|m_{\mathrm{ref},\mathrm{u}}m_{\mathrm{ref},\mathrm{l}}|<15\mathrm{mrad}`$. This ensures that the track segments are not affected by tubes in which delta rays444Energetic primary electrons knocked from the tube walls or gas atoms by the incident charged particle. Delta rays can lead to smaller measured drift times if they are emitted in the direction of the anode wire. distorted the drift time measurement. The remaining events are called good tracks. For the determination of wire positions, tracks with high muon momentum are desirable to limit the effects from multiple scattering. Thus, the estimated muon energy is required to be larger than $`2.5\mathrm{GeV}`$. The efficiencies of these track cuts are given in table 1.

Tubes in the test chamber with hits generated by delta rays would affect the measurement of the wire position, and are therefore rejected: It is demanded that the tube have contributed to a track segment in the test chamber. Also, a loose cut on the difference between the radius prediction from the reconstructed track segments and the measured drift radius is made. For the track segment reconstructed in the test chamber this cut is at 0.7 mm, while it is at 1.5 mm for the predicted radius from the reference track segments extrapolated into the test chamber555The larger value for the reference chamber track segments is due to the long extrapolation distance into the test chamber.. The efficiencies of these cuts are given in table 2. They are very close to 100% and thus do not bias the wire position measurement.

### 8.2 Determination of the Wire Positions

The reconstructed reference chamber track segments are extrapolated into the test chamber and the drift radius measured in the test chamber is compared with the prediction from those track segments. The weighted average $`r_{\mathrm{ref}}`$ of the drift radius predictions from the two track segments

$$r_{\mathrm{ref}}=\frac{(1/\sigma _\mathrm{u}^2)r_{\mathrm{ref},\mathrm{u}}+(1/\sigma _\mathrm{l}^2)r_{\mathrm{ref},\mathrm{l}}}{(1/\sigma _\mathrm{u}^2)+(1/\sigma _\mathrm{l}^2)}$$
(9)
is taken, where each uncertainty $`\sigma _i`$ ($`i=\mathrm{u},\mathrm{l}`$ for the track segment in the upper and lower reference chamber) depends on the extrapolation distance from the reference chamber and the estimated muon momentum (to account for the uncertainty from multiple scattering). The value of $`\sigma _i`$ is about $`200\mu `$m for the reference chamber closer to the tube hit and ranges from 350 to 800 $`\mu `$m for the farther chamber.

If a given wire is displaced by $`\delta y`$ in the $`y`$ direction and $`\delta z`$ in the $`z`$ direction, the measured drift radius $`r_{\mathrm{drift}}`$ in the tube will show a systematic deviation $`\mathrm{\Delta }r`$ from the predicted reference radius,

$$\mathrm{\Delta }r=r_{\mathrm{drift}}r_{\mathrm{ref}}=\delta ym\delta z+𝒪(m^2),$$
(10)
as shown in figure 5.

### 8.3 Determination of the Wire Position Perpendicular to the Chamber Plane

The displacement $`\delta z`$ of an individual anode wire is determined from the slope in a linear fit to the distribution of $`\mathrm{\Delta }r`$ vs. $`m`$.

The $`z`$ displacement of an entire tube layer and its tilt corresponding to a rotation about the $`x`$ axis are determined from the distribution of measured $`\delta z`$ values of all tubes as a function of their $`y`$ coordinate666Here, $`y`$ is taken to be the nominal wire position. Because of the large chamber width, this approximation does not bias the result.. A linear fit to this $`\delta z`$ distribution is performed separately for each layer. The offset and slope of the fitted line are taken to be the $`z`$ displacement and tilt of the layer, respectively.

In the further analysis, each wire is assigned the measured $`z`$ position of the layer at the nominal $`y`$ wire position. The single wire position measurement is used for wires where it deviates from the layer position by more than 3 standard deviations.

### 8.4 Determination of the Wire Position in the Chamber Plane

Using the $`z`$ position measurement described above, the displacement $`\delta y`$ of an anode wire is obtained as

$$\mathrm{\Delta }r^{}=r_{\mathrm{ref}}^{}r_{\mathrm{drift}}=\delta ym\delta z^{}\delta y,$$
(11)
where $`r_{\mathrm{ref}}^{}`$ denotes the predicted radius corrected for the wire $`z`$ position and the residual displacement $`\delta z^{}`$ after correction is of the order of the resolution of the $`z`$ position measurement. The average of the track slopes $`\left|m_{\mathrm{all}\mathrm{hits}}\right|_{\mathrm{all}\mathrm{tubes}}`$ is777The deviation from 0 is caused by a slight asymmetry of the trigger efficiency in the $`y`$ direction during data taking. $`0.06`$, the largest value of $`\left|m_{\mathrm{all}\mathrm{hits}}\right|`$ for a single tube being $`0.15`$. Therefore, the systematic error introduced in the last approximation of equation (11) is typically much smaller than $`15\%`$ times the resolution on the $`z`$ position, which does not limit the $`\delta y`$ measurement.

In the chamber plane, the nominal wire positions can be described by a grid of the form

$$y(n)=y_0+ng,$$
(12)
where $`n`$ denotes the index of the tube in the layer, $`y_0`$ the offset of the layer position relative to the other layers in the same chamber, and $`g`$ the mean distance between two neighbouring wires. Fitting the function in equation (12) to the distribution of the wire positions of one tube layer yields the offset $`y_0`$ and the mean distance $`g`$ for that layer.

### 8.5 Test of the Method

In order to test the method described above, one of the few chambers which have been scanned by the X-ray tomograph at CERN has been used as a test chamber. Two tomograph scans of the chamber were performed at a distance of 30 cm from each chamber end. Therefore the positions of the anode wires888Only typically $`80\%`$ of the wires are visible in the tomograph scans because of the support structure of the chamber. are known with a precision of 2 $`\mu `$m.

The analysis uses 28 hours of data taken with this chamber. During this time 1.7 million events per meter along the anode wires ($`x`$ direction) have been recorded. The number of events after each step of the event selection and the number of tube hits in the test chamber used for the wire position measurements are given in tables 1 and 2.

All measurements at the Cosmic Ray Facility are in very good agreement with the tomograph results.

In figure 6, the $`\delta z`$ measurements for single wires performed at the Cosmic Ray Facility are compared with the results of the tomograph scans. The distributions of the difference between the cosmic ray measurement and the tomograph scan have widths of 23 $`\mu `$m and 27 $`\mu `$m. These widths are dominated by the resolution of the Cosmic Ray Measurement Facility.

A better precision can be obtained for the $`z`$ displacement of entire tube layers. A comparison of this measurement and the values for the tilt angle of the tube layers around the $`x`$ axis (the parameters which are more difficult to control during the MDT chamber assembly) with the tomograph results is shown in figures 7 and 8. Here, we achieve a precision of $`4.4\mu \mathrm{m}`$ for the layer shift and $`17\times 10^6`$ for the tilt angle.

The comparison of the $`\delta y`$ measurement with cosmic rays with the tomograph results shows that the difference between the two measurements has an RMS of $`8.5\mu \mathrm{m}`$ at either end of the chamber (see figure 9). The precision of the Cosmic Ray Measurement Facility is therefore 8.3 $`\mu `$m (after subtraction of the uncertainty on the X-ray tomograph measurement).

The measurements of the layer offset $`y_0`$ and the wire grid constant $`g`$ are compared in figures 10 and 11. A precision of $`1.8\mu \mathrm{m}`$ for the layer offset and $`0.15\mu \mathrm{m}`$ for the grid constant is achieved.

## 9 Conclusions

In the Cosmic Ray Measurement Facility, ATLAS MDT muon chambers are tested and calibrated. The Cosmic Ray Measurement Facility is capable of measuring the positions of the individual anode wires of a chamber with a precision of $`8.3\mu \mathrm{m}`$ in the chamber plane and $`27\mu \mathrm{m}`$ in the direction perpendicular to that plane after 28 hours of data taking.

In addition to the wire positions, the offsets of the tube layers relative to each other in the chamber plane and their grid constants have been determined with precisions of $`1.8\mu \mathrm{m}`$ and $`0.15\mu \mathrm{m}`$, respectively. The positions of the tube layers in the direction perpendicular to the chamber plane and the angles for a rotation around an axis parallel to the wires have been measured with precisions of $`4.4\mu \mathrm{m}`$ and $`17\times 10^6`$, respectively.

At a rate of up to two tested chambers per week, all 88 BOS chambers which are produced in Munich will be tested until the planned start of muon chamber installation in ATLAS.

## 10 Acknowledgements

We would like to thank our colleagues at Dubna and MPI Munich for the excellent collaboration, the X-ray tomograph group at CERN for providing their results, and the staff at LMU for their support. We are pleased to acknowledge the support of the Maier-Leibnitz-Laboratorium of LMU and TU Munich, and the Bundesministerium für Bildung und Forschung, Germany.


format: text

text:
# Vector opinion dynamics in a model for social influence

## 1 Introduction

Much attention has been paid in recent years to the potential applications of the methods of statistical physics to complex systems and phenomena that were traditionally considered to lie far outside the interests of physicists . Modeling such systems as sets of interacting dynamical elements has proven to be a fruitful procedure to capture the essential mechanisms that lead to generic forms of emergent collective behavior –such as pattern formation and synchronization – actually observed in many natural processes. This program has been successfully applied to a large variety of systems in the fields of physics, chemistry, and biology , and has also been extended to the study of social and economical systems . In this class of systems, where elementary agents are human beings, models aim at describing rational behavior, such as decision making or the adoption of a strategy.

In this work we study a generic model of social influence inspired in Axelrod’s proposal , in which the interaction between individuals tends to homogenize the state of the population, in the presence of some pre-existing homogeneity. In Axelrod’s model, each agent is characterized by a “cultural” state defined as a string of *features*, each one of them being a *trait* of integer value. The agents are arranged on a bidimensional lattice, and can interact with nearest neighbors. The dynamical rules of the model prescribe that two agents can interact if they have some cultural similarity, v.g. if some of their features have equal traits. This pre-existing similarity defines a probability that allows the active agent to further approach the state of her partner, by adopting one of her neighbor’s cultural features. The system may evolve into an inhomogeneous state of “cultural regions” whose number and sizes depend on the parameters of the model. Interesting analysis are given in Refs. , where statistical properties of the dynamics are shown.

In this context, we use the term opinion for the set of individual attributes that are subject to social influence and can be updated as a result of the interaction between agents. Opinions can be modeled, in a simple approach, as binary numbers  and their dynamics have been first studied in economy . Axelrod’s “cultural” state, restricted to features made up of binary traits (0 or 1), corresponds to what we call here “opinion” state, and has been used in Refs. . It may be interpreted as the answers to a survey with dichotomic (yes/no) questions aimed at defining the agent’s position on a given subject. Alternatively, opinions can be considered as continuous variables. The conditions to reach consensus in this case were analyzed in several works .

In our model, we use opinions defined by a binary string, i.e., the state of each agent is a vector of components 0 or 1. Since we are interested in the propagation of opinions within a population, agents are randomly selected in pairs, instead of the regular lattice of Axelrod’s model, which is more appropriate for the geographical spread of culture. If the opinions of the selected agents are close enough, they interact. The probability of interaction is the same for all agents, and is controlled by a threshold-like parameter. This is also in contrast with Axelrod’s model, where the probability of interaction increases strictly with the cultural similarity. As for the underlying structure of the population, we analyze two cases. The first corresponds to the situation of complete mixing, where any pair of individuals can interact and adjust their opinions. In the second, we assume that the agents are in the nodes of a small-world network and the interaction process occurs only between connected agents.

In the following section we introduce the model for the opinion dynamics for the complete mixing case, and study the dynamics of the system. In Sec. III, the properties of the stationary state are analyzed. In Sec. IV, we introduce the small-world model and its main results, and emphasize similarities and differences between both cases. Finally, we summarize and discuss our most important results.

## 2 Social influence and the dynamics of opinion

Our model consist of a set of agents that update their opinions as a result of mutual influence. We consider a population of $`N`$ agents, each of them characterized by a binary vector of $`k`$ components representing her opinion. The number of possible different opinions is $`2^k`$, and $`k`$ is the dimension of the opinion space. The initial condition is a random and uniform distribution of opinions. At each time step $`t`$ two randomly chosen agents meet. If the difference between their opinions is lower than a threshold $`u`$, one of them copies each opinion component from the other with probability $`\mu `$. This event is called an “interaction” in the rest of the paper. Clearly, its occurrence is governed by $`u`$, the parameter that characterizes the existence of some –even if partial– previous coincidence of opinions. The probability of change of an agent’s opinion is $`p=1(1\mu )^d`$, where $`du`$ is the distance between the opinions of the interacting pair. Thus, $`\mu `$ plays the role of a convergence parameter of the model.

To quantify the difference between two opinions we use the Hamming distance, i. e., the number of different components between the two binary vectors. In consequence, the interaction process occurs when agents agree in at least $`ku`$ opinion components. The threshold $`u`$ is taken as constant in time and across the whole population.

In our simulations, we consider $`16N4096`$ and $`k=10`$, so that $`2^k=1024`$. Systems with $`N<2^k`$ are said to be “diluted,” whereas systems with $`N>2^k`$ are “saturated,” since there will always be at least two agents with the same opinion. Our numerical calculation thus cover both extremes of saturation. We have used several values of the convergence parameter, between $`\mu =0.066`$ and $`\mu =1`$. We show here the results corresponding to these two values, that serve well as characterization of the two observed dynamics. As we show below, the number of time steps needed to reach the stationary state, $`t_f`$, strongly depends on $`\mu `$ and scales with $`N^2`$. We take $`t_f=500N^2`$ for the lowest value $`\mu `$ used and $`t_f=100N^2`$ for the highest value.

Our analysis begins with the study of opinion evolution from a random distribution in a well mixed population, where any two agents are allowed to interact. We find that the system evolves, at long times, toward clusters of homogeneous opinions. As a characterization of the state of the system we compute, at each time step, the number of agents with the same opinion and calculate the number of clusters of different opinions, $`N_c`$. In Fig. 1 we plot the evolution of $`N_c`$ for several values of the threshold. Figure 1(a) shows the results for the convergence parameter $`\mu =0.066`$, and Fig. 1(b) those for $`\mu =1`$. A low value of $`\mu `$ insures a much slower approach to a stationary state. At short times we find that $`N_c1`$, as expected for a system with a random distribution of opinions. At long times, for sufficiently large thresholds, the number of clusters is $`N_c=1`$, i.e. full consensus is achieved. If the threshold is low enough we have, instead, $`N_c>1`$ at long times. This means that clusters of different opinions exist in the stationary state, in spite of the fact that complete mixing of the interaction mechanism and the homogeneity of the initial configuration allow the propagation of an opinion through the whole population. The transition from a multi-clustered state to a single cluster one occurs at higher values of $`u`$ for the fast $`\mu =1`$ system than for the slower $`\mu =0.066`$ one. This transition is also sharp for $`\mu =0.066`$ (Fig. 1(a)) and smooth for $`\mu =1`$ (Fig. 1(b)), a point to which we will return below.

Figure 1(a) shows, moreover, that the evolution in the number of clusters is not monotonous. There are randomly distributed events where $`N_c`$ increases, caused by the opinion rearrangement that takes place before two clusters collapse. This feature disappears when averaging over different realizations of the system, as is the case of the rest of the results shown in this paper.

Another characterization of the opinion dynamics is provided by the evolution of the cumulative number of interactions in the system, $`N_i(t)`$. If the threshold is large (of the order of $`k`$), every pair of agents selected in a given time step does interact, i.e. $`N_i(t)=t`$. In the opposite limit, if $`u`$ is zero there are no interactions and $`N_i(t)=0`$ for all $`t`$. For intermediate values of $`u`$, the quantity $`N_i(t)`$ illustrates the dynamical state of the system. In Figs. 2(a) and 3 we plot $`N_i`$ as a function of time for different values of the threshold $`u`$, and for $`\mu =0.066`$ and $`\mu =1`$ respectively.

In Fig. 2(a) we find three regimes. An initial one is characterized by a linear growth of $`N_i`$, with a slope that depends on the threshold $`u`$. This regime is relatively short, and can be seen more clearly in the double logarithmic plot in the inset of Fig. 2(a), where the vertical separation of the lines represent the difference in their slopes. This transient regime corresponds to the first stage of clusterization. Since the initial condition is uniform over the phase space, every agent finds her randomly picked potential partner at a random distance, uniformly distributed over the phase space. Due to the fact that the interaction becomes effective only if the partner lies within a ball of radius $`u`$ from the agent, the probability of interactions is proportional to the fraction of the phase space covered by such balls. If $`uk`$, these balls cover almost all the space and the activity is the fastest possible, with $`N_it`$, as can be seen in the figure, for the case $`u=9`$. For smaller $`u`$, the activity is slower and the slope of $`N_i(t)`$ tends to 0. This initial regime continues during the first stage of the clusterization process, and Fig. 2(b) allows a visualization of its extent. Here we show the complement of the number of clusters, $`2^kN_c(t)`$, for $`k=10`$, corresponding to the simulation run of $`u=3`$ shown in Fig. 2(a). A small arrow indicates the end of this regime, where the clusterization reaches a plateau. At this point, the system is composed of a great number of clusters distributed over the phase space. Since the convergence parameter $`\mu =0.066`$ allows an interacting pair to approach their opinions without collapsing them into a single one, most of these clusters lie close to each other. This is in contrast to what we find for $`\mu =1`$, as discussed below. After this plateau, the clusters have become near enough to interact more often, and an acceleration of the activity appears for all values of $`u`$. This is apparent in Fig. 2(a) as an upward bending of $`N_i(t)`$, and in Fig. 2(b) as a nearly linear growth of $`2^kN_c`$. This mechanism rapidly reduces the number of clusters, until there are only a few of them remaining (typically around 10 for $`u=3`$). This marks the beginning of the third and final regime, indicated by the second arrow for the case $`u=3`$ in both figures. The remaining clusters are very close to each other, and the dynamics proceeds at the maximum possible speed of one interaction per time step until the final state of a single opinion is reached, for all $`u>2`$.

In Fig. 3 we show the curves $`N_i`$ vs. $`t`$ that correspond to the highest convergence parameter $`\mu =1`$. For all values of $`u`$ we find two regimes, a linear one at short times with a slope dependent of $`u`$ (for the same underlying reasons as in the case $`\mu =0.066`$), and a long time behavior with a slower increase with time. After the initial linear regime in which there is an intense activity, the interactions become more spaced in time and finally stop, when the stationary state is reached. The reason for this slowing down, at variance with the acceleration displayed when $`\mu =0.066`$, resides in the extreme value of the convergence parameter. Since $`\mu =1`$, the opinions of a pair of interacting agents become identical after the interaction. The clusters, in consequence, become effectively isolated from one another, and intermediate opinions do not exist after the first clusterization stage . Eventually, and depending on the value of $`u`$, a final state with several clusters becomes stationary.

In Fig. 4 we plot the average over $`100`$ realizations of the total number of interactions $`N_i`$ needed to achieve the stationary state, as a function of the system size $`N`$ for two values of the threshold $`u`$, $`k=10`$ and $`\mu =1`$. We observe an approximately linear increase of $`N_i`$ with $`N^2`$, with a slope depending on the threshold, as shown in the inset. For all system sizes we find the same dependence with $`N`$, which coincides with the one observed for the time $`t_f`$ needed to reach the stationary state. This means that the ratio $`N_i/t_f`$, measuring the average number of interactions per time unit during the complete evolution of the system, is independent of the system size. For $`\mu =0.066`$ the dependence on $`N`$ is the same, even though there is a factor between 10 and 100 in the total number of interactions needed to reach the final state. Since the behavior of $`N_i`$ does not depend on $`N`$ in a qualitative way, we use an intermediate representative value, $`N=2^k`$, in the rest of the paper.

## 3 Characterization of stationary states

As shown in the previous section, a stationary state is achieved after a number of interactions that depends on the convergence parameter $`\mu `$ and on the particular realization. We are interested now in the dependence of the final state of the system on the threshold and the convergence parameter.

In Fig. 5 we plot the number of clusters in the stationary state, $`N_c`$, as a function of the threshold $`u`$ for two values of the convergence parameter, $`\mu =0.066`$ and $`\mu =1`$. Each dot of the curve corresponds to an average over $`100`$ realizations. In both cases the value of $`N_c`$ for $`u<2`$ is very similar and corresponds to a situation where several clusters of different opinions coexist in equilibrium. In the region $`u>7`$, both cases are also similar, with $`N_c=1`$. This is a situation where a single opinion always prevails. In the region $`2<u<7`$ the mean number of clusters depends on the convergence parameter. The transition from the multi-cluster state to the single opinion state is smooth for $`\mu =1`$, whereas for $`\mu =0.066`$ there is a sharp transition at $`u3`$. For the system size shown here this last transition is found for $`\mu <0.99`$. A signature of the deep difference between the two transitions is their behavior in a finite size analysis. We show this in Fig. 6 for $`\mu =0.066`$ and $`\mu =1`$, for a range of system sizes that correspond to increasing saturation, since $`k`$ remains fixed. It is clear that the transition becomes sharper with increasing system size when $`\mu <1`$, and that it preserves its smoothness when $`\mu =1`$. This is an indication that there is a further transition at $`\mu =1`$ between the two types of behavior. Only the fastest systems ($`\mu =1`$) allow the persistence of multi-opinion populations at intermediate values of the threshold $`u`$. In Figs. 6(a) and (b) we show a diversity of behaviors that correspond to the different regimes of saturation. This analysis does not allow a further characterization of the transitions. Indeed, the thermodynamic limit would require not only that $`N\mathrm{}`$ but also that $`k\mathrm{}`$, at a fixed saturation in opinion space. Such systems, with opinions defined by infinite features, go beyond the interest of the present study. In a related context, this has been carried out in Ref. , where the switch from the monocultural to the multicultural states in Axelrod’s model is identified as a phase transition. The roles of noise and of network structure have been analyzed in Refs. .

The transition as a function of the threshold can be further characterized analyzing the distribution of the population of clusters, $`P(x)`$, for different values of $`u`$. If the system is in a homogeneous state, all the agents having the same opinion, then $`N_c=1`$ and the distribution has a peak at $`xN`$. A situation of many clusters with a few agents in each cluster is characterized by a distribution $`P(x)`$ with a dominant peak in $`x1`$. In Fig. 7 we show the results for the case $`\mu =0.066`$ and two values of $`u`$ near the transition. For $`u<2`$ the histogram has a single peak in $`x=1`$ whereas for $`u>3`$ it has a single peak in $`x=N`$. The threshold at which the transition between a clusterized state and a homogeneous state takes place is between $`u=2`$ and $`3`$. This result coincides with that of Fig. 5.

A very different sequence of histograms is obtained for $`\mu =1`$, as we show in Fig. 8. The transition from a highly clustered state (the histogram with a peak in $`x=1`$, at $`u=2`$) to a homogeneous state (the histogram with a peak in $`x=N`$, at $`u=7`$), takes place gradually for intermediate values of $`u`$. This result is also in agreement with the behavior observed in the $`N_c`$ curve.

Note that to describe completely the state of the system we need both quantities, $`N_c`$ and the sequence of histograms. As an example, observe that in Fig. 5, $`N_c4`$ for both $`u=2`$ with $`\mu =0.066`$ and $`u=4`$ with $`\mu =1`$. However, the histograms are very different in Fig. 8 with $`u=2`$ and Fig. 9 with $`u=4`$.

The results of this section indicate that a well-defined transition between a homogeneous state and a highly clustered state takes place in this model as a function of $`u`$. A complementary result is that the parameter $`\mu `$ not only modifies the convergence speed towards equilibrium but also changes the structure of the stationary state. This last result contrasts with Ref. , where the authors report that the convergence parameter only modifies the time needed to reach the equilibrium.

## 4 Complex mixing vs. small-world network

The two previous sections have dealt with the case of complete mixing, in which any agent can interact with any other one in the system. How do the dynamical and steady states change if the underlying structure of the population is different? Here we analyze the case in which the agents are situated at the nodes of a small-world network , which represents more accurately some aspects of real social populations. At variance with a complete mixing system, a small-world network exhibits some degree of local clustering in the neighborhood of each agent. We introduce the parameter $`p`$, the reconnection probability, that measures the randomness of the small-world and interpolates between an ordered lattice ($`p=0`$) and a random graph ($`p=1`$). Another property that characterizes the small-world is the network connectivity $`K`$. The small-world network is built from a one-dimensional regular lattice in which each node is linked to its $`2K`$ nearest neighbors. Then, each link is rewired with probability $`p`$ to a randomly chosen node. Multiple links are forbidden and disconnected networks are discarded. Every agent in the small-world network has, on the average, $`2K`$ nearest neighbors.

The interaction process occurs only between linked sites. Now, we randomly choose an agent and then we select, also at random, one of his neighbors to interact with. The mechanism of the interaction is the same as in the complete mixing model. To ease the comparison with previous results, we use $`N=1024`$ and $`k=10`$ in all the figures of this section.

Figure 9 shows the average number of clusters $`N_c`$ as a function of the threshold $`u`$ for a small-world network with $`K=3`$, $`\mu =1`$ and three different values of the reconnection probability $`p`$. We compare this results with those obtained in the complete mixing case of Fig. 5 (open circles). The behavior of $`N_c`$ for $`p0`$ is very similar to the results of the complete mixing case in all the range of thresholds. The only qualitatively different situation is the case $`p=0`$, that corresponds to the agents in a regular lattice. In this case, and for $`u<7`$, the number of clusters with different opinions is significantly larger than those corresponding to other values of $`p`$. The regular lattice –and only the one which is completely regular– is less effective at disseminating the opinions. The same dependence with $`p`$ was observed for the case of $`\mu =0.066`$ (not shown here).

These results indicate that the dynamics in a small-world configuration does not differ from that of a well mixed system, and the only different case is the one corresponding to $`p=0`$, actually, not a small-world network but a regular lattice. The transient evolution of single realizations of such a system is shown in Fig. 10, where we plot $`N_c(t)`$ for three values of the threshold. The parameters are the same as in Fig. 1, where we use the lowest value of $`\mu `$ to have the time evolution in detail. Note that the system on the lattice needs roughly $`10`$ times more iterations than the completely mixed one. We also observe the random fluctuations originated by opinion rearrangements in the clustering process.

Finally, in order to analyze the effect of the connectivity in the system, we change the values of $`K`$. As can be observed in Fig. 11, the mean number of clusters tends to the complete mixing case at high values of $`K`$. The number of nearest neighbors necessary to recover that behavior depends on the threshold.

## 5 Conclusion

We have studied a generic model of social influence in which agents adjust their opinions as a result of random binary encounters whenever their difference in opinion is below a given threshold.

We have analyzed the dynamical behavior of the system and found that the cases $`\mu =1`$ and $`\mu <1`$ behave differently. For the extreme case $`\mu =1`$ we observe an initial linear regime of intense activity that ceases when intermediate opinions disappear. This first stage is followed by a different regime in which the interactions become more spaced in time because the opinion clusters are effectively isolated, even for $`u`$ very near $`k`$. Finally, a saturation occurs when the stationary state is reached. The slower case $`\mu =0.066`$ shows three regimes. The first stage of clusterization is similar to the one observed in the previous case, but it produces a great number of clusters that lie close to each other. A second stage starts, in which the clusters interact more often, giving rise to an acceleration of the activity until a few clusters remain in the system. Finally, in the third stage the dynamics proceeds at the maximum speed due to the interaction between clusters that are close enough. The dynamical process stops when a single opinion is reached (for $`u>2`$) or several isolated clusters are formed (for $`u<2`$). We can see, then, that the speed of the convergent dynamics plays a relevant role not only in the final composition of the population but also in the intermediate stages of evolution.

We also found that both the total number of interactions between individuals and total time needed to reach the equilibrium, increase quadratically with the size of the system. Consequently, the ratio between these two quantities, $`N_i/t_f`$, is independent of the system size. We obtained this result for all the values of $`\mu `$ used.

As for the stationary state properties, we have found that the interaction between individuals restricted by a proximity threshold results into clustering of opinions with a number of clusters decreasing with the increase of the threshold. Our model presents a transition driven by the threshold, from a state with homogeneous opinion (consensus) to a phase in which several clusters with different opinions coexist. The nature of the transition depends on the convergence parameter $`\mu `$, that changes the structure of the stationary state. This was also observed in the sequence of histograms of population of clusters near the transition. Finite size scaling supports the existence of such transition.

Finally, we have studied the influence of the underlying structure on the behavior of the system. We found that small-world networks and the complete mixing case have a very similar behavior, except for the limiting case of regular graphs ($`p=0`$ in the small-world model). This last case is less effective to disseminate the opinions than the others. The peculiarity of the case $`p=0`$ is found also in other contexts, within the study of small-world networks. There has been found that the small-world phenomenon appears at any value of the network disorder $`p`$, for sufficiently large systems. Geometrical properties such as the average path length , and some dynamical ones, such as the ferromagnetic properties of an Ising model , have been shown to behave critically at $`p=0`$. This is not the case in some other dynamical models, such as the propagation of an epidemic , and of a rumor , that behave critically at a finite value of $`p`$ in the thermodynamical limit. In our model of opinion dynamics, the effect of disorder is even stronger, since any nonzero value of $`p`$ induces a behavior that is indistinguishable from both the fully disordered $`p=1`$ case and the complete mixing case. We also found that greater network connectivity gives rise to a decrease in the number of clusters. The number of connected neighbors needed to reach the complete mixing ones depends on the threshold.

Several open problems associated with opinion dynamics are worth mentioning. For instance, the interplay between long-range effects, such as mass-media broadcasting of individual or collective opinions, and local interactions is an essential mechanism in modern opinion formation. Boundary effects and inhomogeneity conditions, including spatial variation of population density, may also strongly affect the dynamics of social influence and the properties of the resulting stationary state. Finally, the effects of noise, communication errors, and in general random fluctuations, should not be disregarded in any realistic representation of information spreading. The present model is expected provide a versatile tool to analyze these questions.

## Acknowledgement

M.F.L. thanks the Solid State Group of Centro Atómico Bariloche for the computational facilities.


text:
# Character of Atomic Vibrations in a Lennard-Jones Glass

## I Introduction

Vibrational properties of disordered media have been reviewed by various authors . Subsequent to most of these reviews, Grest, Nagel, and Rahman (GNR)  did a careful characterization of the normal modes of the Lennard-Jones (LJ) glass. Our previous numerical studies, motivated by heat conductivity experiments and by theoretical controversy, focussed on vibrations in amorphous Si (Am-Si) . Here we present results for a very large LJ glass model, with local atomic structure very different from Am-Si. These results confirm the GNR picture. We use the term “glass” simply to mean an amorphous solid, and the abbreviation am-LJ to denote the LJ glass.

The model studied is a monatomic $`6`$-$`12`$ Lennard-Jones system, $`V(r)=4ϵ[(\sigma /r)^{12}(\sigma /r)^6]`$, with potential parameters suitable for Argon: $`ϵ/k_B=125.2`$ K and $`\sigma =0.3405`$ nm. This potential gives a satisfactory description of crystalline Ar . In its face-centered cubic (fcc) structure, each atom is surrounded by 12 nearest neighbors located a little inside the minimum of the $`V(r)`$ at $`r=2^{1/6}\sigma `$. The second shell is more distant by $`2^{1/2}`$.

Standard molecular dynamics simulations were done on systems with $`N`$=$`500`$, $`2048`$ and $`6980`$ atoms, enclosed in a cubic box with periodic boundary conditions. A truncated ($`r_c=2.6\sigma `$) LJ potential was used. The simulated density was $`\rho =42\times 10^3`$ mol/cm3, corresponding to the reduced value $`\rho ^{}`$=$`\rho \sigma ^31`$. After an equilibration run at $`T=400`$ K using NVE molecular dynamics, gradual cooling was performed down to a temperature ($`T=150`$ K) slightly above the melting temperature. Finally, conjugate gradient minimization found the glassy minimum configuration.

The resulting structure of am-LJ is somewhat like that of a smeared-out argon crystal. The density is about 5% smaller. There may be icosahedral aspects in the local packing , but our models are more similar to the “LJ glass” of Simdyankin et al.  than to their “IC glass” with strongly icosahedral features. A shell of radius r contains on average $`(4\pi /3)\rho r^3`$ atoms, where $`\rho `$ is the number density. The mean distance from an atom to its $`N`$’th neighbor is thus $`(3N/4\pi \rho )^{1/3}`$.

For am-LJ we calculated the mean distance to the $`N`$’th neighbor, and the variance $`<\delta r^2>^{1/2}`$ of this quantity, which are plotted in Fig. (1). The closest 10 neighbors cluster in a narrow interval near the distance 3.7Å found in crystalline Ar. The 11th to 13th neighbors have less well-defined positions. The 15th to 18th atoms have the least well-defined positions, lying near the second neighbor shell of the crystal. More distant neighbors start to be packed as in a random gas, with distance starting to scale like $`N^{1/3}`$.

Normal modes of vibration were calculated in harmonic approximation, using the same periodic boundary conditions used to compute the positions. The calculations were done on Njal, a Beowulf cluster with 81 dual pentium III processors running at 0.933 GHz. The Scalapack routine PDSYEV was used. The largest system size, 6,980 atoms, corresponds to a real-symmetric dynamical matrix of dimension 20,940, the largest that fits conveniently into Njal’s memory. This matrix required about 100 minutes to diagonalize. The vibrational density of states is shown in Fig. (2). All three models give essentially the same result, which is shown here for the 6980-atom model, and was first found by Rahman et al. . The results for the glass are compared with the density of states of the crystal, which closely resembles experiment on crystalline Ar . However, the LJ glass has a surprisingly different density of states, being featureless, with no distinction between longitudinal and transverse modes.

To pursue this issue, it is interesting to look at the “bond-stretching parameter” $`S_a`$ of the $`a`$’th normal mode, defined as

$$s_a=\left[\frac{_{}|(\stackrel{}{u}_i^a\stackrel{}{u}_j^a)\widehat{n}_{ij}|^2}{_{}|(\stackrel{}{u}_i^a\stackrel{}{u}_j^a)|^2}\right]^{\frac{1}{2}}$$
(1)
where $`\stackrel{}{u}_i^a`$ is the displacement eigenvector on the atom with position label $`i`$ and position $`\stackrel{}{R}_i`$. The unit vector $`\widehat{n}_{ij}`$ points in the direction $`\stackrel{}{R}_i\stackrel{}{R}_j`$. In other words, $`S_a`$ measures for mode $`a`$ the relative amount of stretching of each near-neighbor bond. For the crystal, this is rigorously defined by using the 12 bonds from each atom to its first neighbors. For the glass, we sum all neighbors whose separation is less than the mean distance to the 12.5th neighbor. The results are plotted in Fig. (3).

This figure shows a very striking difference between the crystal and the glass. Crystalline normal modes can be labelled by the wavevector $`\stackrel{}{k}`$. In Ar, there are only three normal modes at each $`\stackrel{}{k}`$. They have usually quite different bond-stretching characters, and can normally be classified as longitudinal (L) or transverse (T). This leads to a highly differentiated spectrum of states, seen in the figure as a broad range of values of $`S_a`$ at each frequency. In the glass, the bond-stretching character $`S_a`$ collapses to a narrow interval whose width diminishes to zero as $`N`$ increases to infinity. Lacking the quantum number $`\stackrel{}{k}`$, each mode has totally mixed character. This property was noticed previously for other glasses . In an infinite glass, with infinitely many modes at each frequency, one could make a unitary rotation (within the subspace of fixed frequency) to a basis where certain modes would look approximately pure L in some local region. However, at long enough distances this character is completely lost. At low frequencies, purified modes defined this way could be chosen have a fairly sharply defined local value of $`\stackrel{}{k}`$. At large enough distances, the mode would evolve into a spherical mixture of different directions of $`\stackrel{}{k}`$. As frequency increases, there is a crossover (the “Ioffe-Regel” (IR) crossover) to the regime where the local $`\stackrel{}{k}`$ can hardly be defined because it decays in a distance of one or two neighbors. The IR crossover and other low-frequency properties will be the subject of a future study.

It was noticed by GNR that a given normal mode, insofar as it was L, expressed its L character at specific wavevectors, and insofar as it was T, expressed that character at other wavevectors. This enabled them to draw pseudo-dispersion curves showing very different character for L and T aspects. This does not contradict our view that a given mode cannot be classified. Each mode is simultaneously L and T, probably even at the same point in space unless purified as described above.

In principle, normal modes of vibration, just like electron eigenstates in quantum theory, can be rigorously classified as extended (E) or localized (L). In three dimensions the vibrational spectrum has sharp E/L boundary frequencies (“mobility edges,” $`\omega _c`$) separating these two kinds of modes. The IR crossover is low in the spectrum of a glass, but the E/L boundary is usually high in the spectrum.

An hypothesis that the mobility edge coincides with the IR crossover  has not proved true for any model that has been studied numerically, but has remained a popular misconception. Our large system provides a useful confirmation that previous studies have correctly located the mobility edge. We compute the participation ratio $`p_a`$ for each normal mode, defined as

$$p_a=\frac{[_i|\stackrel{}{u}_i^a|^2]^2}{_i|\stackrel{}{u}_i^a|^4}.$$
(2)
The values of $`p_a`$ are plotted on a logarithmic scale versus eigenfrequency for all three models in Fig. (4).

Notice first that the value of $`p_a`$ hovers at slightly higher than 0.5$`N`$ for a significant part of the spectrum (frequencies 2.5 to 6.5 meV). These are very well delocalized states whose amplitude does not vary too much from atom to atom. At low frequencies (below 1 meV for the 6980-atom model, and pushing up to 2 meV for the 500-atom model) there is a trend to lower values of $`p_a`$, and there occur some idiosyncratic modes with $`p_a`$ reduced by 10 or more. Their interpretation remains difficult . At least some of this behavior has to be attributed to finite-size artifacts. We hope to make a detailed study later. Finally, at high frequencies, very low values of $`p_a`$ occur. These states are strongly Anderson localized. The mobility edge is estimated to be at $`\omega _c=8.2\pm 0.1`$ meV. This estimate was found by fitting a straight line to the logarithmic data at the point where $`p_a0.25N`$. Extrapolating this line back to the horizontal line $`p_a=0.5N`$, the intersection occurs at around 8.2 meV, with uncertainty which diminishes as the system size increases. This value of $`\omega _c`$ corresponds to having the lower 93% of the spectrum delocalized, and the upper 7% localized, and agrees with the result of GNR .

These results make an interesting comparison to our earlier work on amorphous Si . One major difference is that the density of vibrational states in am-Si retains much of the structure seen in the crystal, with acoustic and optical types of vibrations separated by a region of low density of states. To the contrary, in am-LJ, the density of states loses differentiation and is completely featureless. Apart from this, other properties of am-LJ vibrations are reminiscent of those in am-Si. We see similar behavior in the collapse of the bond-stretching character onto a single frequency-dependent curve. The anomalous states at very low frequency, and the strongly localized states at very high frequency, have similar character. In both glasses, fewer than 10% of the modes are localized, always at the upper end of the spectrum.

## ACKNOWLEDGMENTS

We thank G. S. Grest, J. Fabian, J. L. Feldman and G. Viliani for help. Work at Stony Brook was supported by NSF grant DMR-0089492. We thank K. K. Likharev for the use of the Njal computer cluster, which was purchased using a grant from DOD’s DURIP program.


text:
# Effects of intermediate bound states in dynamic force spectroscopy

## I Introduction

The last decades have witnessed a revolution in the methods to observe and manipulate single bio-macromolecules or bio-molecules complexes. New micromanipulation techniques have especially been put forward to probe the folded structure of proteins and to quantify the strength of adhesion complexes titin ; Poi01 ; Sim99 ; Nis00 ; Pie96 . An important step in this direction is the proposal of the group of Evans to use soft structures to pull on adhesion complexes or molecules at various loading rates (dynamic force spectroscopy or DFS) Eva . Moving the other end of the soft structure at constant velocity induces on the complex a pulling force that increases linearly in time $`f=rt`$. Measuring the typical unbinding time $`t^{}`$ yields an unbinding force $`f^{}=rt^{}`$ that depends on the pulling rate $`r`$. The typical outcome of such experiments is a plot of $`f^{}`$ vs. $`\mathrm{log}r`$ composed by a succession of straight lines with increasing slopes (force spectrum). It has been argued and that it is possible to deduce the value of some relevant structural parameters of the complex by analysing force spectra thanks to an adiabatic Kramers picture. This picture consist in considering the unbinding process as the thermally activated escape from bound states over a succession of barriers along a one-dimensional path crossing a mountainous energy landscape Eva . Within this scheme, each straight line of the force spectrum witnesses the overcome of an energy barrier and its slope maps the barrier to a distance $`x`$ along the pulling direction. This procedure has been shown to yield reasonable values for a few systems, and has been conforted by numerical simulations Hey00 .

Subsequently, theoretical studies have refined the above original model, e.g., by inclusion of rebinding events Sei02 , study of time dependent loading rates Eva , or consideration of more complex topographies Str00  and topologies of the energetical landscape Bar02 .

In the present paper we explore the potential influence of the existence of intermediate bound states on the experimental dynamic response of adhesion complexes as probed in DFS. To achieve this goal, we first revisit the analysis of the escape from a bound state consisting of an arbitrary number of barriers along a 1D path under the application of an external load (in line with earlier studies of Strunz et al. Str00 ), and then discuss the implications of this analysis for the interpretation of experimental data. In (III) the standard picture is recalled, together with its two underlying assumptions. In (IV), we first relax the a priori assumption of a deep fundamental bound state and provide a general expression that relates the typical rupture force to the loading rate (within a single escape rate approximation). The practical implications of this new formula \[Eq. (15)\] are discussed, and in particular we comment upon intrinsic ambiguities in inferring informnation from a $`[\mathrm{log}(r),f^{}]`$ plot. Then, we show in (V) that in the presence of multiple bound states it may be necessary to relax the other assumption (a single typical rupture force for each loading rate) as multimodal rupture force distributions naturally show up, a feature recently observed in lipid extraction experiments Eva02 .

## II Model and Notations

Figure 1 illustrates the energy landscape of a one-dimensional escape path with $`N`$ energy barriers and wells. The position and energy of the $`i`$th energy well ($`i=0`$ marking the fundamental bound state, and $`1iN1`$ the intermediate ones) are denoted by $`x_i`$ and $`E_i`$, respectively. Similarly, the position and energy of the $`j`$th energy barrier are denoted by $`\widehat{x}_j`$ and $`\widehat{E}_j`$, respectively (where $`1jN`$). For convenience, without loosing generality, we set $`x_0=0`$ and $`E_0=0`$ for the fundamental bound state. The unbound “state” is on the right hand side of the $`N`$th barrier. If the energy differences $`(\widehat{E}_iE_i)`$ and $`(\widehat{E}_{i+1}E_i)`$ exceed $`k_\mathrm{B}T`$ the transition rates $`k_i^{}`$ (and $`k_i^+`$) from the $`i`$th energy well over the left $`i`$th barrier (and right $`(i+1)`$th barrier, respectively) can be written according to the Kramers formula

$`k_i^{}`$ $`=`$ $`\omega _0\alpha _i\widehat{\alpha }_i\mathrm{e}^{(\widehat{E}_iE_i)/k_\mathrm{B}T},`$ (1)
$`k_i^+`$ $`=`$ $`\omega _0\alpha _i\widehat{\alpha }_{i+1}\mathrm{e}^{(\widehat{E}_{i+1}E_i)/k_\mathrm{B}T},`$ (2)
where $`\omega _0`$ is a typical attempt frequency, $`\alpha _i`$ and $`\widehat{\alpha }_j`$ are geometric factors characterizing the shape of the $`i`$th energy well and $`j`$th energy barrier, respectively. Note that there is no transition from the fundamental bound state to the left, therefore, $`k_0^{}0`$.

We assume throughout the paper that the energy wells and barriers are sharp, so that for any loading force $`f`$ their locations remain constant, and their energies change as $`E_i(f)=E_i(0)fx_i`$ and $`\widehat{E}_j(f)=\widehat{E}_j(0)f\widehat{x}_j`$. To simplify the notations, wherever the argument of the energies and transition rates is omitted, a loading force $`f`$ is implicitly assumed.

Finally it will prove convenient to introduce a few compact notations. For any $`0i<jN`$ we denote the distance between the $`i`$th well and the $`j`$th barrier (on the right) by $`\mathrm{\Delta }x_{i,j}=\widehat{x}_jx_i`$, and their energy difference by $`\mathrm{\Delta }E_{i,j}=\widehat{E}_jE_i`$. We also define a formal (effective) rate constant from the $`i`$th well over the $`j`$th barrier on the right as

$$k_{i,j}=\omega _0\alpha _i\widehat{\alpha }_j\mathrm{e}^{\mathrm{\Delta }E_{i,j}/k_\mathrm{B}T}.$$
(3)
Obviously $`\mathrm{\Delta }E_{i,j}`$ and $`k_{i,j}`$ implicitly depend on $`f`$, whereas $`\mathrm{\Delta }x_{i,j}`$ are constants given the assumption of the previous paragraph.

## III Stochastic kinetics of unbinding under external forces: standard description and corresponding approximations

We first recall the standard description of the ”force spectrum”, which relies on two major assumptions, namely the single escape rate (SER) and the deeply bound fundamental state (DBFS) approximations.

### III.1 Single escape rate approximation

If for any experimentally relevant load $`f`$ the equilibration within the bound states is much faster than the escape to the unbound state, the unbinding can be described by a single load dependent escape rate $`k(f)`$. Following the calculation of Evans Eva01 , if rebinding is negligible (which is the case in most experimental situations), the probability $`P(t)`$ of remaining in the bound state at time $`t`$ (the survival probability of the bond) then decreases as

$$\frac{\mathrm{d}P(t)}{\mathrm{d}t}=k(rt)P(t).$$
(4)
The solution of this differential equation is $`P(t)=\mathrm{exp}[_0^tk(rt^{})dt^{}]`$. The probability density for unbinding between times $`t`$ and $`t+\mathrm{\Delta }t`$ is $`p_t(t)=\mathrm{d}P(t)/\mathrm{d}t=k(rt)P(t)`$, from which, after changing the variable from $`t`$ to $`f`$, one gets the probability density for the distribution of the unbinding force: $`p_f(f)=(1/r)k(f)P(f/r)`$. The typical unbinding force $`f^{}`$ is defined as the peak of this probability density: $`\mathrm{d}p_f(f)/\mathrm{d}f|_{f=f^{}}=0`$, which yields the simple formula

$$\frac{\mathrm{d}\tau (f)}{\mathrm{d}f}|_{f=f^{}}=\frac{1}{r},$$
(5)
where $`\tau (f)1/k(f)`$ denotes the load dependent mean escape time. This formula gives the loading rate $`r`$ at which the typical unbinding force is $`f^{}`$. For practical purposes it is often necessary to invert this relation to, e.g., predict the typical unbinding force for an experimentally imposed loading rate.

To set a reference for further comparison, we explicitely invert the above relation in case of a single barrier, i.e. when $`\tau (f)=\tau (0)\mathrm{exp}(f\widehat{x}_1)`$, and obtain

$$f^{}=\frac{k_\mathrm{B}T}{\widehat{x}_1}\mathrm{ln}\left[\frac{r\tau (0)\widehat{x}_1}{k_\mathrm{B}T}\right].$$
(6)
As mentioned in the introduction, the escape over a single barrier results in a single straight line in the force spectrum. The experimental observation of a linear segment consequently gives hints as to the structure of the energy landscape, in particular the slope of the segment permits to deduce a distance $`\widehat{x}_1`$ between the energy well and the barrier.

### III.2 Deeply bound fundamental state approximation

Assuming further that the fundamental bound state is much deeper than the intermediate ones: $`E_i(f)E_0(f)k_\mathrm{B}T`$ for any experimentally relevant load $`f`$ (i.e. before unbinding has statistically almost certainly occurred, see Fig. 2 a.), Evans has shown that the mean escape time from the fundamental bound state to the unbound state is well approximated by Eva01 ; Eva02

$$\tau (f)=\underset{j=1}{\overset{N}{}}\frac{1}{k_{0,j}(f)}=\underset{j=1}{\overset{N}{}}\frac{\mathrm{e}^{f\mathrm{\Delta }x_{0,j}/k_\mathrm{B}T}}{k_{0,j}(0)}$$
(7)
This allows one to obtain an explicit $`r`$ vs. $`f^{}`$ relationship by plugging Eq. (7) into Eq. (5), which yields the compact formula

$$r=\left[\underset{j=1}{\overset{N}{}}\frac{\mathrm{\Delta }x_{0,j}}{k_\mathrm{B}T}\frac{\mathrm{e}^{f^{}\mathrm{\Delta }x_{0,j}/k_\mathrm{B}T}}{k_{0,j}(0)}\right]^1.$$
(8)
This equation predicts a spectrum $`f^{}`$ vs. $`\mathrm{log}(r)`$ consisting of a succession of at most $`N`$ segments of increasing slopes, each of which yielding an information $`\mathrm{\Delta }x_{0,j}`$ about an intermediate barrier.

## IV Beyond the deeply bound fundamental state approximation

In this section we relax the approximation made in the last subsection (IIIB), generalize accordingly equations (7) and (9), and discuss the experimental implications of this generalization.

### IV.1 Refined theory

In general, it is possible that for large enough forces one or more of the intermediate bound states become deeper than the fundamental bound state before unbinding has occured (see Fig. 2 b.). In such cases the above DBFS approximation breaks down. However, we show below that it is still possible to compute rather simply the escape time $`\tau `$ from the ”bound state” to the ”unbound state”, provided we maintain the assumption of a single escape rate $`1/\tau (f)`$.

Let us put the system into its fundamental bound state, and let it evolve according to the transition rates given in Eqs. (1) and (2). Whenever the system gets into the unbound state (by making a transition over the outermost barrier) let us place it back into the fundamental bound state. The stationary state of an ensemble of such systems is characterized by a probability current, which is constant everywhere and equal to $`1/\tau `$ by definition. To calculate $`\tau `$ we have to solve the following system of equations:

$`P_ik_i^+P_{i+1}k_{i+1}^{}`$ $`=`$ $`1/\tau 0iN2,`$ (9)
$`P_{N1}k_{N1}^+`$ $`=`$ $`1/\tau ,`$ (10)
$`{\displaystyle \underset{i=0}{\overset{N1}{}}}P_i`$ $`=`$ $`1,`$ (11)
where $`P_i`$ denotes the probability of being in the $`i`$th bound state ($`0iN1`$). The first $`N`$ equations describe the probability current over each of the $`N`$ barriers, and the last equation is just the normalization condition. These $`N+1`$ linear equations uniquely determine the $`N+1`$ variables ($`P_i`$ and $`\tau `$), and can be solved easily in a recursive way. First, $`P_{N1}\tau `$ can be expressed from Eq. (10), and then $`P_{N2}\tau `$, $`\mathrm{}`$, $`P_0\tau `$ recursively from Eq. (9) yielding

$`P_i\tau `$ $`=`$ $`{\displaystyle \frac{1}{k_i^+}}+{\displaystyle \frac{k_{i+1}^{}}{k_i^+k_{i+1}^+}}+\mathrm{}+{\displaystyle \frac{k_{i+1}^{}\mathrm{}k_{N1}^{}}{k_i^+k_{i+1}^+\mathrm{}k_{N1}^+}}`$ (12)
$`=`$ $`{\displaystyle \underset{j=i+1}{\overset{N}{}}}{\displaystyle \frac{1}{k_{i,j}}},`$
where Eqs. (1), (2), and the definition (3) have been used. Note that because the $`k_{i,j}`$ are only formal definitions, constructed as products and ratios of the single-barrier rates (1) and (2), they are meaningful even if $`\mathrm{\Delta }E_{i,j}<0`$. From the normalization (11) one can easily express $`\tau `$ as

$$\tau =\underset{i=0}{\overset{N1}{}}P_i\tau =\underset{i=0}{\overset{N1}{}}\underset{j=i+1}{\overset{N}{}}\frac{1}{k_{i,j}}.$$
(13)
The sum is dominated by the smallest effective rates, which are the bottlenecks of the unbinding process. Consequently, this formula remains a good approximation for $`\tau `$ even if some of the barriers disappear at big loads, because the corresponding formal transition rates make negligible contributions. By indicating the load force $`f`$ explicitly, we arrive at

$$\tau (f)=\underset{i=0}{\overset{N1}{}}\underset{j=i+1}{\overset{N}{}}\frac{1}{k_{i,j}(f)}=\underset{i=0}{\overset{N1}{}}\underset{j=i+1}{\overset{N}{}}\frac{\mathrm{e}^{f\mathrm{\Delta }x_{i,j}/k_\mathrm{B}T}}{k_{i,j}(0)}.$$
(14)
which generalizes (7). An analytic formula can be given for the $`f^{}`$ vs. $`r`$ relationship by plugging Eq. (14) into Eq. (5) :

$$r=\left[\underset{i=0}{\overset{N1}{}}\underset{j=i+1}{\overset{N}{}}\frac{\mathrm{\Delta }x_{i,j}}{k_\mathrm{B}T}\frac{\mathrm{e}^{f^{}\mathrm{\Delta }x_{i,j}/k_\mathrm{B}T}}{k_{i,j}(0)}\right]^1,$$
(15)
This generalization of equation (8) is one of the main results of this paper. Let us briefly comment on immediate features of this new formula.

First, Eq. (8) is easily recovered from (15) assuming a DBFS. Indeed, the assumption $`E_i(f)E_0(f)`$ implies $`k_{0,j}k_{i,j}`$, if $`i>0`$ \[see Eq. (3)\] and therefore, the relation $`P_i/P_01`$, if $`i>0`$ is deduced from Eq. (12). The probability to find the system in the fundamental bound state is close to $`1`$. So, the sum over $`i`$ in Eq. (14) is dominated by the contributions of the effective escape rates from the $`0`$th well only. Finally, the sum over $`i`$ (labeling the intermediate states) is reduced to its sole first term too, and Eq. (15) becomes identical to Eq (8).

Second, each of the $`N(N+1)/2`$ terms of Eq. (15) alone would yield a straight line in the $`f^{}`$ vs. $`\mathrm{log}(r)`$ plot. However, at any loading rate the highest force value (the uppermost segment, corresponding to the most difficult transition) limits the unbinding process, therefore, the $`f^{}(r)`$ curve is expected to closely follow the upper envelope of these segments \[see Fig. 3(a)\]. Depending on the position of the lines, this upper envelope can consist of up to $`N(N+1)/2`$ linear segments.

Third, this last point is clearly at odds with the prediction within the DBFS approximation. Indeed assuming a DBFS corresponds to forbidding the display in the force spectrum of the $`N(N1)/2`$ segments corresponding to the probing of the escape from an intermediate bound state (see Fig. 2. b.).

### IV.2 Practical implications: Ambiguity in the determination of ”structural” parameters

We now insist on some practical implications of the above general description. We do not attempt a full inspection of all the possible dynamic responses of arbitrarily complex systems, but rather focus on two simple examples in order to stress that the main features of the energy landscape can in general not be unambiguously inferred from $`[\mathrm{log}(r),f^{}]`$ plots. To emphasize the experimental relevance of this discussion, we use for the parameters values comparable to those observed in experimental systems. Specifically, we take the geometric factors $`\alpha _i`$s and $`\widehat{\alpha }_j`$s to be all equal to $`1`$, $`\omega _0=10^8`$ s<sup>-1</sup> and $`k_\mathrm{B}T=4\times 10^{21}`$J.

#### IV.2.1 Ambiguity in determining the barriers positions

Fig 3 (a) and Fig 3 (b) display two force spectra as obtained from Eq. (15). Both correspond to energy landscapes with two barriers. Though the two $`[\mathrm{log}(r),f^{}]`$ plots are almost identical they are related to very different sets of values for the energy levels and positions (along the pulling direction) of the wells and the barriers.

Fig 3 (a) corresponds to the situation where the standard picture to account for the two segments is well suited Eva . At low force, the escape from the fundamental $`0`$th state over the outermost barrier is the limiting process. The slope of the first segment is proportional to $`k_BT/\widehat{x}_2`$. For the highest forces (above $`30\mathrm{p}\mathrm{N}`$) the energy of the external barrier is reduced below $`\widehat{E}_1`$ and the deepest bound state remains located at $`x_0=0`$. The process that mostly impedes the unbinding is the overcome of the innermost barrier $`\widehat{E}_1`$ with a rate $`k_{0,1}`$. The slope of the curve is now larger and proportional to $`k_BT/\widehat{x}_1`$.

Fig. 3 (b) corresponds to an energy landscape for which the above explanation is inappropriate. At low force the unbinding kinetic is controlled by the escape from the fundamental state over the outermost barrier again. But, for pulling forces larger than $`30\mathrm{p}\mathrm{N}`$ this outer barrier remains the highest \[see inset in Fig 3.(b)\]. However the slope of the spectrum increases as in the (a) case. The reason is that the deepest (and most occupied) bound state is now located at $`x=x_1`$ and the presence of the second segment actually witnesses the escape from this intermediate state to the unbound state with a rate $`k_{1,2}`$. The value of the second slope scales therefore with $`k_BT/(\widehat{x}_2x_1)`$. Since the escape rate $`k_{1,2}`$ in the (b) case is equal to the escape rate $`k_{0,1}`$ in the (a) case, the two spectra in Fig. 3 turn out to be indistinguishable and cannot be a priori associated with one of the two possible landscapes.

#### IV.2.2 Ambiguity in determining the number of barriers

After having shown with the simple example above that ambiguity can exist in determining distances from dynamic force spectra, we show here that even more strikingly it is impossible in general to assess the number of wells and barriers. Again we use a simple example to do so.

Fig. 4 displays two force spectra obtained using Eq. (15). They are both well approximated by a succession of three segments with increasing slopes. Again, the two $`[\mathrm{log}(r),f^{}]`$ curves are very similar although they are constructed from landscapes that do not even comprise the same number of peaks and wells.

In Fig. 4 (a) the three segments describe the escape from the same fundamental state over the three distinct energy barriers. The larger the pulling force the closer the limiting barrier to the fundamental state \[see inset in Fig. 4 (a)\].

In Fig. 4 (b), the landscape consists of only two barriers. However, the force spectrum reveals that three different escape processes can limit the unbinding kinetic. At low forces ($`f50\mathrm{p}\mathrm{N}`$) the two observed linear segments results from the escape form the fundamental state over the two peaks at $`\widehat{x}_1`$ and $`\widehat{x}_2`$ respectively. Conversely, at high forces it is the escape from the deeply lowered intermediate state over the outer barrier that determines the escape rate (see inset in Fig. 4 (b), drawing with dash-dot line). With the chosen parameters the effective rate $`k_{0,3}`$, $`k_{0,2}`$ and $`k_{0,1}`$ in Fig. 4(a) case correspond respectively to $`k_{0,2}`$, $`k_{0,1}`$ and $`k_{1,2}`$ in Fig. 4(b) case. Thus the two plots are indistinguishable and cannot be used to predict the number of barriers along the 1D escape path.

In conclusion of this subsection, we suggest great care in inferring features of the underlying energy landscape from dynamic force spectroscopy experiment. Our generalized equation may be helpful in dealing with the corresponding ambiguity as it allows (with some work) to generate various landscapes that can account for the observed data, whereas Eq. (8) can only yield a single set of parameters ( e.g. those used for the plots in Fig. 3 (a) and Fig. 4 (a)).

## V Beyond the single escape rate approximation: Multimodal unbinding force distributions

Up to this point we have been considering a generalized theory in which the deeply bound fundamental state (DBFS) approximation is dropped, but the unbinding is still approximated as a simple first-order escape process. Indeed, the validity of Eq. (15) relies on the assumption that at any moment the distribution of the populations of the bound states can be well approximated by the distribution corresponding to a homogeneous stationary current.

This is, however, not always the case. As we stated earlier, the sum of the $`1/k_{i,j}(f)`$ terms in Eq. (14) is dominated by the smallest effective rate constant $`k_{i^{},j^{}}(f)`$ corresponding to the slowest effective transition. A consequence of this is that all the bound states located to the left of barrier $`j^{}`$ are close to equilibrium (because of the slow outflow over barrier $`j^{}`$), and the population of any state located to the right is negligible (because they practically belong to the unbound state). Now, If the slowest transition rate changes from $`k_{i^{},j^{}}(f)`$ to $`k_{i^{\prime \prime },j^{\prime \prime }}(f)`$ as the loading force $`f`$ is increased, and if $`j^{\prime \prime }<j^{}`$, then a considerable population might remain in the intermediate bound states between the new and the old limiting barriers, $`j^{\prime \prime }`$ and $`j^{}`$ respectively. This residual population is incompatible with the new stationary current dominated by $`k_{i^{\prime \prime },j^{\prime \prime }}(f)`$, and must escape in a different way, yielding a secondary maximum of the unbinding force distribution (see Fig. 5b).

The escape of the majority of the population (located to the left of the new limiting barrier $`j^{\prime \prime }`$) can still be characterized by Eq. (14) of our generalized theory. On the other hand, we have to slightly modify this formula to describe the escape of the residual population (trapped between the new and old limiting barriers). Since $`j^{\prime \prime }`$ is the limiting barrier now, almost the entire residual population can escape without ever jumping backward over barrier $`j^{\prime \prime }`$. Therefore, for the residual population we can consider barrier $`j^{\prime \prime }`$ as a reflecting boundary, and describe the escape by our general theory in this modified potential. Eq. (14), e.g., changes accordingly:

$$\tau (f)=\underset{i=j^{\prime \prime }}{\overset{N1}{}}\underset{j=i+1}{\overset{N}{}}\frac{1}{k_{i,j}(f)}=\underset{i=j^{\prime \prime }}{\overset{N1}{}}\underset{j=i+1}{\overset{N}{}}\frac{\mathrm{e}^{f\mathrm{\Delta }x_{i,j}/k_\mathrm{B}T}}{k_{i,j}(0)}.$$
(16)
Consequently, the absolute maximum of the unbinding force distribution always follows the upper envelope of the $`N(N+1)/2`$ lines, however, some secondary maxima might also appear at lower forces, which follow the upper envelope of only a subset of the lines \[comprising $`(Nj^{\prime \prime })(Nj^{\prime \prime }+1)/2`$ elements\]. Such secondary maxima of a multimodal force distribution give important information on the internal structure of the energy landscape of the unbinding path, and makes the determination of the number and positions of the energy wells and barriers less ambiguous. It is actually a nice achievement of our generalized theory to be able to make sense of the segments of secondary maxima in a unique frame for fitting parameters (see e.g. Fig. 5, where the segment corresponding to the secondary maximum corresponds to the transition from the well 1 over the barrier 2, a step neglected in the DBFS approximation). The possibility of a bimodal distribution for the case of a two-state system has already been reported by Strunz et al. Str00 , and our description systematizes and generalizes their findings.

To provide a simple illustration for the somewhat formal discussion above, we also focus on a system consisting of two bound states, as depicted in Fig. 5 (a). Increasing the force, the limiting transition rate changes from $`k_{0,2}`$ to $`k_{1,2}`$ and then to $`k_{0,1}`$. In the range of the loading rate $`r`$ between about $`10^4`$ and $`10^5`$ pN/s the intermediate bound state (1) has enough time to accumulate a large population, which is then flushed by the $`k_{1,2}`$ transition before the transition $`k_{0,1}`$ flushes the rest from the fundamental bound state (0). In the range above $`10^6`$ pN/s the intermediate bound state (1) cannot accumulate much of the population, but it still possesses a small fraction of the initial equilibrium distribution, which is again flushed by the $`k_{1,2}`$ transition first.

Very recently the group of Evans actually reported the experimental occurrence of a bimodal force distribution Eva02 . The corresponding experiment consisted in pulling on ”diC14 PE” lipids from a bilayer made of ”C18:0/1 PC” lipids. With the help of our generalized theory, Evans and Williams were able to fit their data and interpret the results in terms of an energy landscape with two barriers \[personal communication, see Ref.  in Ref. Eva02 \].

## VI Conclusion

In this paper, we have revisited the standard theory used to account for the dynamic response of molecular stickers. Our refined description, valid for an arbitrarily complex one-dimensional energy landscape, has allowed us to highlight several practical consequences of the diversity of the possible unbinding scenarios. For example several markedly different energy landscapes can yield the same rupture force distribution. To resolve this ambiguity other experimental technics, e.g. flow chamber experiments Pie02 , are then required. We have also identified the physical origin of multimodal unbinding force distributions and shown how their analysis provides informations on the unbinding pathways.