Matthias Büchner

Jean-Philippe Loisel

The sonorous beat of optical waves in a Michelson interferometer

This paper describes time dependant phase modulation experiments using an optical Michelson interferometer. Changes in the path length of the interferometer arms induce phase shifts in time and, thanks to a loudspeaker, the interferometric signal is audible. We produce time fluctuating refractive indices and mechanical motions of one reflecting mirror. A sinusoidal motion is applied, thus creating a phase modulation of one interferometer arm wave. Sidebands of the laser frequency appear which interfere with the other interferometer wave. Our hearing sense is quite developed and we use it to analyse the non linear character by the acoustic timbre. Finally we transfer information by phase modulation. Our set-up benefits to the visually-impaired community as temporal phase changes can be analysed by the auditory sense, but this unconventional use of an optical interferometer also attracts a general audience.

I. INTRODUCTION

The Michelson interferometer is a standard instrument to demonstrate interference of an optical wave. It can also be used as a precision instrument to detect tiny changes in the length of optical interferometer paths. One or both interferometer output beams can be sent to photodetectors and the optical intensity signals can be recorded, allowing high precision measurements of laser wavelengths [START_REF] Dobosz | Overview of the laser-wavelength measurement methods[END_REF][START_REF] Pedregosa-Gutierrez | Computer-controlled high-precision Michelson wavemeter[END_REF] and distance changes [START_REF] Riles | Gravitational waves: Sources, detectors and searches[END_REF] . One impressive example is the detection of gravitational waves with the LIGO experiment, rewarded in 2017 by the Nobel prize for R. Weiss, B.C. Barish and K.S. Thorne [START_REF] Weiss | LIGO and the Discovery of Gravitational Waves, I: Nobel Lecture[END_REF][START_REF] Barish | LIGO and Gravitational Waves II: Nobel Lecture[END_REF][START_REF] Thorne | LIGO and Gravitational Waves, III: Nobel Lecture[END_REF] , Traditionally in education, the fringe patterns are observed visually and indeed their beauty is impressive and flawless.

Such an interferometer can also be used to demonstrate interferences in time, i.e. beats of optical waves. Mickey McDonald et al. [START_REF] Mcdonald | Visible optical beats at the hertz level[END_REF] described a Mach-Zehnder interferometer with an acousto-optical modulator (AOM), placed in each interferometer arm. These AOM shift differently the laser frequency by about 80 MHz and their frequency difference is some Hertz. The optical beating of the two laser waves can be directly observed: the intensity at any point within the pattern is seen to oscillate sinusoidally in magnitude. The beauty of this experiment is the observation by eye of a pulsation in the Hertz region, at low frequencies never directly introduced in the experiment.

Our approach also aims to observe beat notes of optical waves, but our set-up places the beat observation into the sound world. In fact, our auditory sense [START_REF] Stix | How Hearing Works [Video[END_REF][START_REF] Rigby | Human brains naturally tuned to hear music[END_REF][START_REF] Hawkins | Human ear[END_REF][START_REF] Hartmann | Signals, Sound, and Sensation[END_REF][START_REF] Plack | The Sense of Hearing[END_REF] can discern precise frequencies: our fascination with music is based on this fact [START_REF] Freke | Sine language[END_REF][START_REF] Hartmann | Principles of Musical Acoustics[END_REF] . Our set-up doesn't need AOMs and their highly stabilised power supplies. The use of the auditory sense is of interest for the visually-impaired community but also adds an amusing educational aspect for a wide-spread audience.

a) Electronic mail: matthias.buchner@irsamc.ups-tlse.fr

We add another component to a standard optical Michelson interferometer: we monitor the optical interference signal by a photodiode and we send its amplified signal to a loudspeaker or headphone. We placed one mirror on a piezo actuator controlled translation stage. Our goal is to detect relative optical path length changes; i.e. phase shifts in time, by simply listening to the sound.

Changing optical path length in time can be created by moving the position of the piezo-mounted mirror. Another possibility is to alter the refractive index of the optical medium probed by one interferometer arm. In our case, the optical medium is air, its pressure and temperature fluctuations change the optical path difference thus the signal amplitude. Interference in time, i e. beating of the optical waves takes place and thanks to our set-up this beat is audible.

In the following, we recall briefly the theoretical background. The experimental set-up will be described. We will then detail some experiments :

1. Touching the mirror mounts. This induces tiny changes in the mechanical path length difference, i.e. phase shifts in time which can be heard.

2. Making audible temporal changes of the refractive index of the air.

3. Applying a sinusoidal phase modulation by modulating the position of one mirror. Harmonics are created which depend on the motion amplitude.

The order and intensity of the harmonics can be analysed by listening to the sound as the acoustical timbre changes.

4. Finally the electrical signal of a MP3 player is amplified and sent to a piezo mounted mirror, acting on the mirror position: The sound is transmitted by phase modulation and the interference signal contains this information.

Finally we will give hints to lower the costs, i.e. to replace some commercial material by low-cost home made components. A conclusion completes our paper.

II. THEORETICAL BACKGROUND

We assume, that the light source emits a plane wave, with frequency ν L and wavelength λ L . The laser wavelength bandwidth is supposed to be negligibly small. The electric field is then given by the well-known equation:

E(t, x) = E 0 cos ϕ(x, t), with (1)
ϕ(x, t) = 2π[ x λ L -ν L t]
The incident wave is split into two parts with equal intensities (see fig. 1). We use a 50%-50% beamsplitter, its coated side points to the right. One wave is travelling to the right; the other is reflected at the second (coated) surface of the beamsplitter and travels upwards. Both beams are reflected by two mirrors and are superposed on the beam splitter, which recombines both beams. The upper wave is phase shifted by 4πL/λ L , where L is the distance between the beam splitter and the mirror. A phase shift φ 1 is added to include other differences in the optical path, e.g. due to perturbation of the air refractive index. The other wave undergoes a similar phase shift of 4π(L+∆x)/λ L +φ 2 , where L+∆x is the distance between the right mirror and the beamsplitter. The wave ("1") travelling downwards passes through the beam splitter, while the other wave ("2") is reflected at the right (air-side) surface, undergoing a π phase shift. We dropped all phase shifts related to the beamsplitter thickness. Both waves are superposed and form an output beam, situated at the lower part where we place our detector. Defining the wave phases ϕ 1(2) at the detector position y D :

ϕ 1 = 2π(y D + 2L)/λ L + φ 1 -2πνt ϕ 2 = 2π(y D + 2L + 2∆x)/λ L + φ 2 + π -2πνt
We can write Ẽ(t) 2 = [E(t, y D )/E 0 ] 2 /2:

Ẽ(t) 2 = [cos ϕ 1 (t) + cos ϕ 2 (t)] 2 (2) = cos 2 ϕ 1 + cos 2 ϕ 2 + 2 cos ϕ 1 cos ϕ 2 = cos 2 ϕ 1 + cos 2 ϕ 2 + cos (ϕ 1 -ϕ 2 ) + cos (ϕ 1 + ϕ 2 )
The detected intensity, I(t), is proportional to ⟨ E(t, y D ) 2 ⟩ T , where < . > T stands for the average over the detector response time T with T >> 2π/ν L . The terms ⟨ cos 2 ϕ 1(2) ⟩ T are oscillating at twice the optical frequency ν L =474 THz (λ L = 632.8 nm), and each term averages to 1/2. For the same reason, the interference term ⟨cos (ϕ 1 + ϕ 2 )⟩ T vanishes. We get A beam splitter separates this wave into two equal parts, one travelling to the right, the other upwards. Thanks to two mirrors (Thorlabs PF10-03-P01), both waves are reflected and superposed on the beam splitter.

I(t) ∝ 1 + cos (ϕ 1 -ϕ 2 ) (3) ∝ 1 + cos (4π∆x/λ L + φ 2 -φ 1 ) + π
The beam splitter recombines the two waves and two output beams can be distinguished, one at the bottom of the figure, the other on the left, with the outgoing beam pointing towards the laser. Remember, that this outgoing beam should not be injected into the laser. Large fluctuations in intensity and frequency or even laser damage may occur otherwise.

We monitor the first output beam using the amplified signal from a photodetector (Thorlabs PDA36A2). The signal is connected to a oscilloscope, a loudspeaker and to a computer via an analog-digital convertor. We used the software Labview (National Instruments) to communicate with a multifunction card (National Instruments: PCI 6259).

The right mirror is mounted on a commercial translation stage (ACMEL TM-21-10). A piezo actuator (Piezomechanik GmbH Pst 150/5/40) with a maximum extension of 40 µm was placed between the moving element and an end-mounted positioning screw.We use a low voltage which is composed by a DC bias and a time varying voltage. This pilot voltage is amplified by a homemade amplifier up to 100 V to drive the piezo actuator. The choice of all components was dictated by the availability of the material in our laboratory.

At the beginning, it is convenient to align the Michelson interferometer as usual, i.e. in the traditional optical way. Using a ruler, we equalize the distances between the mirrors and the beamsplitter. We superposed the two outgoing beams till fringes are visible. Then, we act carefully on the inclinations of the mirrors to get the fringes as big as possible. The observed signal of the photodetector is given by:

I(t) = I 0 [ 1 + V cos ( 4π∆x(t) λ L + φ 2 -φ 1 ) ] (4)
with I 0 the mean value, proportional to the laser power, and ∆x(t) the optical path length difference of the interferometer arms. This difference can be varied by the piezo actuator. V is the fringes visibility with V ≤ 1.

There are many reasons which lead to smaller V, we can cite non-perfect superposition of the interfering beams and a path difference longer than the coherence length.

Once fringes are visible, a triangular signal from a waveform generator drives the position of the mounted mirror in a linear forth-back motion. The photodetector signal is monitored by an oscilloscope. We used the modulation signal to improve the alignment.

Typical values for our Michelson interferometer are V= 95 % and I 0 = 4.5 V. Our maximal phase sensitivity ∆I/∆ϕ = I 0 V is about 4.3 mV/mrad and for the displacement sensitivity we get ∆I/∆x= 85 mV/nm. I(t) contains two parts (see eq.( 4)). The constant part, I 0 , can not be heard. The other part is proportional to the fringe visibility V and I 0 and this part produces the sound. For typical speeds of the mirror mount around v M = 70 µm/s, ν M is about 221 Hz and this interference signal can be heard thanks to the loudspeaker.

IV. APPLICATIONS

In this section, we will discuss several applications of our set-up and detail the underlying physics. The first two and the last examples (subsections A,B,D) concern qualitative experiments, while C is of more quantitative nature.

A. Transforming mechanical motion into sound: playing like a disk-jockey

A "starter" to demonstrate the sensitivity of the interferometer is to transform mechanical vibration into sound. By changing the DC bias voltage, we set our interferometer to work around mid-fringe, i.e. at maximum phase sensitivity. Just by pushing gently on the fixed mirror, one interferometer path length is changed in time, the fringes move and you can create sounds in a disk-jockey way. i.e. moving vinyl discs by hand back and forth.

The file "disk-jockey.wav" 15 contains one example of sound recording. During the first 7 seconds, one can hear the vibrations induced by the footsteps. Then, the fixed mirror is pushed and the time varying path length difference leads to intensity fluctuation, thus to audible sounds.

One might think, that the pushed mirror changed the direction of the laser beam and only a fraction of the laser beam hit the sensitive area of the photodetector thus inducing a different signal. While continuing the pushing, we blocked at 14 seconds the right interferometer arm: no sounds are created. Unblocking this interferometer arm, again the sound is audible. We now blocked the other arm and again the sound died. This experiment proves that both interferometer arms are needed to create sounds. We also pushed the other mirror with the same results (29 seconds).

An alternative approach is to reduce the visibility to zero by misaligning the interferometer. The resulting signal is given only by I 0 and no sounds are audible.

B. Listening to refractive index fluctuations

We can induce refractive index fluctuations on one interferometer arm. The refractive index of air depends on the temperature and the pressure. Taking the equation of Birch and Downs [START_REF] Birch | An Updated Edlén Equation for the Refractive Index of Air[END_REF] for dry air at temperature 20 • C and pressure 1 bar (100 kPa), a pressure change of +10 kPa increases the refractive index of ∆n ≈ 3×10 -5 (λ L = 633 nm). Over a length of ℓ = 1 mm, this change induces a phase shift of ∆φ = 2π∆n(ℓ/λ L ) ≈ 0.3 rad. For a temperature change of +10 • C, we get ∆n ≈ -0.9 × 10 -5 , i.e. ∆φ ≈ -0.1 rad.

We provoked air density fluctuations in the probe zone. We tried to apply a heating resistance under one interferometer arm, the fringes did move but their motion was too slow to be audible. Better results were achieved with a lightened match, although not to be recommendable nearby optical elements ("allumette.wav" 15 ). We employed an air duster with good results. A weak air stream is directed towards one interferometer beam in a vertical direction.

The file "airduster.wav" 15 contains the sound produced (start at 4 seconds). At 11 seconds, one interferometer arm is blocked, released and at 14 seconds the other arm is blocked. Of course, both interferometer beams should be present to produce the sound. At 19 seconds the air stream is directed towards the blind side of the mirror mount. No sound can be heard, i.e. the sound heard before is not related to a hypothetical mirror vibration produced by the air stream.

C. Generation of frequency harmonics

Our set-up can be used to create frequency combs of the main laser frequency. Theodor W. Hänsch [START_REF] Hänsch | Nobel Lecture: Passion for precision[END_REF] and John Hall [START_REF] Hall | Nobel Lecture: Defining and measuring optical frequencies[END_REF] received the Nobel prize in 2005 for their ground breaking work on laser frequency combs and their important applications in high resolution spectroscopy.

We create frequency combs by phase-modulating one optical wave in the interferometer. A sinusoidal pilot voltage is applied to the piezo actuator, i.e. the mirror is oscillating with an amplitude a m and frequency ν m : ∆x = ∆x 0 + a m sin (ω m t), with the angular frequency ω m = 2πν m . By changing the DC bias voltage, ∆x 0 can be varied slowly in time. Inserting the fast sinusoidal motion, eq. ( 4) reads:

I(t) = I L [1 + V (cos (φ 0 + φ m sin (ω m t)))] (5) = I L [ 1 + V cos φ 0 ( J 0 (φ m ) + 2 ∞ ∑ k=1 J 2k (φ m ) cos (2kω m t)) ) -V sin φ 0 ( ∞ ∑ k=1 2J 2k-1 (φ m ) sin ((2k -1)ω m t) )]
with φ 0 = 4π∆x 0 /λ L , φ m = 4πa m /λ L . We developed I(t) using first-order Bessel functions J k (see e.g. ref. [START_REF] Abranowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] ).

Harmonics of the mirror excitation frequency ν m are created with frequencies and amplitudes governed by J k . If we choose φ 0 = -π/2, the total phase is modulated around zero, i.e. at midfringe. Only odd harmonics exist. For φ 0 = 0, modulation occurs at the maximum intensity. Only even harmonics exist.

We recorded the interference signal of a phase modulation with a frequency ω m = 2π× 220 Hz and a small modulation amplitude of φ m = 0.14 rad. Figure 2 shows our data. We varied slowly the offset of the DC bias voltage, thus scanning the interference fringes. We pass from midfringe through the minimum (17 s recording time) and the maximum intensity (34 s). For this small phase modulation amplitude only the fundamental (k = 1) and the first harmonic (k = 2) contribute to the signal: J 1 (0.14) = 0.06983 (220 Hz) and J 2 = 0.00245 (440Hz). The rectangles on the graph enclose the data subsets taken for two Fourier analysis shown on Figure 2. The upper FFT graph corresponds to data in the mid-fringe window (≈ 39 s) and the lower to data around minimum (≈ 58 s). The first FFT analysis shows the presence of the fundamental frequency only (220Hz) while the second exhibits the second harmonic (440 Hz) only. As expected, these results are in complete agreement with the Bessel function amplitude. We converted the data to a sound file ("220Hz10mV.wav" 15 ).

While scanning the pilot offset voltage, the change in timbre can easily be recognized. At mid-fringe (φ 0 = -π/2, 28s recording time) the fundamental at 220 Hz (music notation: A 3 ) dominates as expected, while for the minima and maxima the second harmonic (440 Hz, A 4 , 17, 34 and 47 s), can clearly be heard. The 440 Hz sound component has lower volume, about 4 % of the 220 Hz one, according to the amplitude ratio |J 2 (0.14)/J 1 (0.14)|.

For the following experiment we set the offset of the pilot function to mid-fringe, i.e. φ 0 = π/2. We applied a 220 Hz modulation changing the modulation amplitude pilot voltage in order to vary the phase modulation amplitude φ m from 0.25 to 7.5 rad in steps of 0.25. The file "midfringes3.wav" contains the recorded sound.

At the beginning the volume of the fundamental frequency increases with higher modulation amplitudes. While increasing φ m , the change in timbre can clearly be recognized. In fact the harmonics get stronger. At 24 seconds the third harmonic (660 Hz, A 5 ) can be distinguished, but this depends on the person listening.

However, one has to be aware that φ 0 drifts slowly as we do not employ any active stabilisation techniques. Additionally, the response of a piezo actuator is non linear to the applied voltage and this is particularly true for large voltage amplitudes. As a consequence, the phase modulation is only approximately linear.

At 49 seconds of recording time, the 220 Hz sound is weak. This corresponds to the first zero of the Bessel function J 1 (φ m ) with φ m = 3.8 rad. Using the first part of eq.( 5), we fitted the interference signal and got φ m = 3.8 rad, which is in agreement with the theoretical value. For high phase amplitude modulations, higher order harmonics are produced, we can speak of a frequency comb.

This frequency comb has frequency spacing in the 100-2000 Hertz region, well below the typical bandwidth of a multimode 633 nm free running He-Ne laser, typically in the GHz region. As a consequence, sidebands can not be resolved in the energy (frequency) spectrum, but as we superpose this frequency comb to a coherent reference wave (the other interferometer arm), frequency fluctuations of the laser source are rejected and beats are produced which are in our case audible. This is the base of homodyne detection technique.

D. Transfer of music by phase modulation of optical waves

An interesting experiment is the transmission of music by phase modulation of optical waves. We adjust our interferometer to mid fringe φ m = π/2. We used a sound file with music ("04-petit-pantin-au-coeur-deglace-eponyme.wav" [START_REF] Laei | Petit pantin au coeur de glace[END_REF] ). The sound amplitudes are converted by a Labview program and a digital-analog convertor produces the driving voltage for the piezo actuator, i.e. sound is converted to mechanical motion. The amplification and the sound volume should be weak in order to privilege the fundamental frequency in phase modulation. Higher amplitudes promote higher harmonics, thus distortion of the sound. The sound file "musicPM.wav" contains the interference signal.

During a demonstration, one person of the audience may be invited to use his/her mp3 player to plug the driving voltage input into the sound port. The sound is transmitted via phase modulation and this can be verified by blocking the interferometer arms one after the other: once again the need of two laser beams for the transmission is demonstrated.

A discussion may be triggered about the advantages and inconveniences of the optical phase modulation techniques. Obviously this transmission protocol is not very handy for every-day use as this setup is very sensitive to vibrations and two coherent laser beams are needed (the two interferometer arms). However phase-modulation, in a radio-frequency or optical domain, is very important for research, we can cite the Pound-Drever-Hall technique for laser locking [START_REF] Drever | Laser phase and frequency stabilization using an optical resonator[END_REF][START_REF] Black | An introduction to PoundDreverHall laser frequency stabilization[END_REF] .

According to quantum mechanics, phase modulation, especially beating (or the transmission of music), is not reserved to acoustic or electromagnetic waves: matter waves can also be phase modulated in time. We can cite the work of Décamps et al. [START_REF] Décamps | Observation of Atom-Wave Beats Using a Kerr Modulator for Atom Waves[END_REF][START_REF] Décamps | Phase modulation of atom waves: theory and experiment using the atom optics analogue of the Kerr effect[END_REF] (and refs. within). The group studied beating of lithium waves by phase modulation in an atom Mach-Zehnder interferometer and transmission of music has been realised by phase modulated atom waves (see Décamps et al. [START_REF] Décamps | Observation of Atom-Wave Beats Using a Kerr Modulator for Atom Waves[END_REF] ).

E. Reducing costs and improvements

The He-Ne laser can be replaced by a low cost laser pointer. The replacement of the electronics and/or the mechanical parts depends on the skills in electronics and mechanics of the user: the high voltage amplification unit and the photodetector may be replaced by home-made versions as well as the mirror mountings. In fact, the pilot voltage is less than 10 V and it must be biased by a static voltage of about 50 V. The piezo actuator can be replaced by an positioner based on electromagnetic forces. Yin-Yen Kuo et al. [START_REF] Kuo | 3D-Printed Linear Positioner with Micrometer Accuracy[END_REF] describe a 3D printed positioner which has micrometer accuracy. The displacement of a permanent magnet is controlled by a coil and the magnet acts on a moving stage.

V. CONCLUSION

We used an optical Michelson interferometer to make time-varying optical path length changes audible. Our hearing, especially recognizing frequency changes, is very sensitive and we use it to demonstrate beat notes of optical waves. We apply this interferometer to recognize mechanical vibration of the interferometer mirrors and refractive index fluctuations of the air.

We applied a sinusoidal phase modulation by mechanically displacing one interferometer mirror. The nonlinear response in phase of the interferometer leads to the creation of harmonics of the fundamental frequency. Optical frequency combs can be created, their beat has timbres as a function of applied phase modulation amplitude. Finally we realized a transmission of music by phase modulation of an optical wave. This non-conventional set-up adds educational aspects, which are original and amusing. Making optical interference audible paves the way for a visually-impaired community to understand interferences of optical waves and but we think that the audible character of interferences also interests a general audience. The "disk-jockey" experiment has already been shown during a public event in our university. The audience comprised pupils of secondary education level. The experiments and the underlying physics described in this paper will serve for educational events, courses at school and talks for general public.

Figure 1

 1 Figure1shows our set-up. The optical part is the well know Michelson interferometer. An initial wave is emitted by the He-Ne laser (Lasos LGR7621, 5mW, λ L =632.8 nm).

FIG. 1 .

 1 photodetector

FIG. 2 .

 2 FIG.2. A sinusoidal phase modulation with φm= 0.14 and ωm = 2π× 220 Hz was applied on one interferometer arm. We varied slowly φ0 as a function of time. The left hand inset shows a zoom on the recorded data demonstrating the applied modulation. A subset of data (inside the rectangles at 39 and 58 s) was selected for a Fast Fourier transform (FFT) analysis. The results are shown in the right hand insets.

ACKNOWLEDGMENTS

We gratefully acknowledge Laurent Polizzi, Eric Panader and Michel Gianesin for the realisation of some mechanical supports. The authors owe particular gratitude to Bruno Lepetit and Stephen Gibson for fruitful discussions and careful review of the manuscript. Financial support from the INP-CNRS and the Université Toulouse III Paul Sabatier is acknowledged.


Benoit Coasne

Annie Grosman

Camille Ortega

M Simon

Adsorption in noninterconnected pores open at one

Keywords: numbers: 68.43.-h, 61.43.Gt, 64.70.Fx, 05.70.Np

come

Adsorption isotherms in mesoporous materials exhibit two main typical features: (i) a sharp increase in the amount of gas adsorbed at a pressure below the saturated vapor pressure of the bulk gas which is attributed to capillary condensation in pores of confining sizes; (ii) the desorption process fails to retrace the path of the adsorption one, so that a hysteresis loop appears revealing the irreversibility of the phenomenon.

Zsigmondy 1 proposed the first explanation of the capillary condensation based on the macroscopic Kelvin equation. In order to explain the hysteresis phenomenon, Cohan 2 went further and proposed different scenari for filling and emptying of a pore which he also described with the macroscopic concept of meniscus. In a cylindrical pore open at both ends, the shape of gas/adsorbate interface is different during the adsorption and desorption processes leading to hysteresis. On the contrary, in a cylindrical pore closed at one end, the adsorption is reversible since the same meniscus is present during the filling and the emptying of the pore.

The irreversibility of the adsorption process in a cylindrical pore open at both ends, expected by Cohan on the basis of thermodynamic considerations, has been predicted by new numerical approaches. Grand Canonical Monte-Carlo (GCMC) simulations 3 , analytical approaches 4 together with Density Functional Theory (DFT) 5 have also allowed proposing that the hysteresis phenomenon is an intrinsic property of the phase transition in a single idealized pore and arises from the existence of metastable gas-like and/or liquid-like states. No calculations have yet been performed for cylindrical pores closed at one end.

Macroscopic descriptions based on Cohan's model are currently used to describe the adsorption in non interconnected mesoporous media and in particular to determine the Pore Size Distribution (PSD) via the Barrett-Joyner-Halenda (BJH) method 6 .

However, the experimental validation of Cohan's model has never been checked since the synthesis of ordered mesoporous media with non interconnected pores open at one or at both ends was unsuccesful for a long time. We report in this paper, for the first time, on an experimental investigation of Cohan's predictions, using mesoporous silicon which exhibits two essential properties: (i) the tubular pores are non interconnected and do not exhibit narrow sections and (ii) it exists under two forms with pores closed at one end when the porous layer is supported by the substrate or open at both ends when the layer is removed from the substrate.

Porous silicon is prepared by electrochemical etching of Si single crystal in a HF solution. Its morphology depends on the type and concentration of the dopant atoms [START_REF] Grosman | Structural and optical prop-erties of porous silicon nanostructures[END_REF] . The adsorption studies presented here have been performed with porous layers prepared from highly boron doped [100] Si substrate. The porosity is controlled by the current density and HF concentration, and the thickness by the duration of the anodic dissolution. The pore volume is determined with high precision by weighing the sample before and after the electrochemical dissolution. To calibrate the porosity and thickness of the layer we need a third weighing after the dissolution of the porous part in a NaOH solution. Note that the porous layers are perfectly reproduced. The porous layers under study has a porosity of 51% and the length of the pores is 21 µm.

The morphology of these layers (see Fig. 1(a) and 1(b)) can be described as honeycomb-like structure with pores perpendicular to the Si substrate separated from each other by Si single crystal walls (inset Fig. 1(b)) of constant thickness (∼ 5 nm).

The polygonal-shaped pores are of the same length and closed at one end by the Si substrate. At the end of the porous layer formation which takes about 20 minutes, an increase in the current density leads to the electropolishing regime during which the Si walls at the bottom of the pores are dissolved so that the porous layer comes off the Si substrate. Such a process takes only a few seconds so that the PSD of the porous layer is not changed. Playing with this possibility, a porous Si layer or a porous Si membrane can be obtained with pores open at one or at both ends, respectively. We have shown that there is no lateral interconnection between the pores in these layers. The method involves deposition by evaporation at ambient temperature, of an aluminum layer, 500 nm thick, on a part (A c and B c ) of the top of two porous samples A and B, followed by a thermal annealing under neutral atmosphere of argon at 450 o C for half an hour to improve the contact between the Al deposit and the silicon. We have shown by Rutherford BackScattering analysis (RBS) with 2 MeV α particles that the Al deposited on the samples is well localized on the top of the pores and constitutes a cap. One of the two samples (A) was thermally oxidized in O 2 enriched in 18 O (99%), at 300 o C, 12 mbar, for one hour. This leads to the formation of a thin silicon oxide (∼ 1 nm) on the silicon walls. The 18 O contents in the porous layers under the Al cap (A c and B c ) and beside it (A o ) were determined by the nuclear reaction 18 O(p,α) 15 N at a proton energy near the narrow resonance at E p = 629 keV (Γ = 2 keV). For the measurements under the Al caps the proton beam energy was increased to take into account the energy loss of the incident protons in the Al layer. The nuclear reaction spectra corresponding to A o and A c are shown in Fig. 2. The α peak from A o indicates the presence in the porous layer of 18 O atoms coming from the thermal oxidation. Due to the presence of the Al cap, the α peak from A c is expected at an energy lower by 100 keV than that from the α peak from A o . The counts from A c and B c which are one order of magnitude lower than that from A o are equal, and hence correspond only to a background coming from nuclear reactions with silicon, boron atoms, etc. This experiment reveals that this medium is composed of non interconnected single pores, the tops of which are in direct contact with the gas reservoir during adsorption. To our knowledge, this is the only porous material for which the absence of interconnection has been experimentally proven.

Another important property of this medium is the constancy of the Si wall thickness revealed by the TEM plane view performed on a porous layer thinned down to 1000 Å, i.e. about 10 times the size of the pores. This is a strong argument against the presence of narrow sections inside the pores, the presence of which would lead, through shadow effects, to an inhomogeneous apparent thickness of the Si walls.

The PSD of the porous Si layers has been estimated from a numerical treatment of the TEM images. The negative TEM photographs are digitized by means of a CCD camera. The value of the threshold to make the image binary is chosen to reproduce the porosity of the sample. Analysis of this binary image allows one to obtain the cross section area and perimeter of each pore. The PSD, which corresponds to cylindrical pores having the same cross section area as the polygonal pores, is shown in the inset of Fig. 1(a). The mean pore diameter is 13 nm with a standard deviation of ±6 nm.

Adsorption measurements were made by a volumetric technique using a commercial apparatus under secondary vacuum of around 10 -6 torr. Nitrogen adsorption isotherms were measured up to the saturating vapor pressure of bulk nitrogen (P 0 ) and at 77 K. Porous Si layers were introduced in a cell just after their preparation. The samples were outgassed at ambient temperature for a few hours. This outgassing does not change the nearly complete passivation of the inner surface by a monolayer of SiH 1,2,3 groups which takes place during the electrochemical etching and which is stable for a few days after the formation [START_REF] Grosman | Structural and optical prop-erties of porous silicon nanostructures[END_REF] .

Figure 3(a) shows N 2 adsorption isotherms at 77 K corresponding to porous Si layers with pores open at one or at both ends, the PSD of which is represented in the inset of Fig. 1(a). These two curves exhibit wide and asymmetrical hysteresis loops with shapes corresponding to type H2 (1985 IUPAC classification) or type E [START_REF] De Boer | The structure and properties of porous materials[END_REF] .

The adsorbed amount of gas sharply increases between 0.55 and 0.8 P 0 and then reaches a plateau region where all the pores are filled with the dense phase. The steep desorption process occurs over a pressure range about three times narrower than that of the adsorption.

Assuming that the density of the dense phase equals that of bulk liquid nitrogen, the pore volume extracted from the value of the adsorbed amount on the plateau region is in very good agreement with that measured by gravimetry.

Cohan's model predicts that, for a cylindrical pore open at both ends, condensation takes place at a pressure P c given by an equation similar to the classical Kelvin equation where the curvature of the meniscus is replaced by a cylindrical curvature :

ln P c P 0 = - σV L RT(r -t(P c )) (1)
where V L is the molar liquid volume, σ is the liquid-gas surface tension at temperature T, r is the radius of the pore and t(P) the thickness of the adsorbed layer at P. For a cylindrical pore closed at one end, surface thermodynamics considerations predict a preferential adsorption on the bottom of the pore; a hemispherical meniscus is then present provided that there is perfect wetting, so that the condensation occurs at a pressure P h given by the following modified Kelvin equation:

ln P h P 0 = - 2σV L RT(r -t(P h )) (2)
The evaporation of the wetting fluid proceeds through the same hemispherical meniscus in both cases i.e. at the same pressure P h . Thus, according to Cohan's model : (i) the adsorption should be reversible in a pore open at one end but irreversible in a pore open at both ends; (ii) the adsorption branch for a pore open at both ends should be located at a pressure (P c ) higher than that for a pore closed at one end (P h ); (iii) for a given size distribution of cylindrical pores open at both ends, the adsorption branch should be steeper than that of the desorption branch (provided that the diameter exceeds 3-4 nm) and the shape of the hysteresis loop becomes rather symmetrical as the pore diameter increases, as shown by numerical calculations based on Eqs. ( 1) and ( 2) in which we have introduced a currently used t(P) proposed by Halsey 11,[START_REF] Gregg | Adsorption, Surface area and Porosity[END_REF] or the t(P) reproducing our experimental values prior to the capillary condensation.

To extract the latter t(P) values, we prepared a porous Si layer with a PSD having a mean value as high as possible (about 50 nm) in order to shift the capillary condensation pressures to the highest possible values. The obtained isotherm exhibits hysteresis loop located above 0.8 P 0 . We have thus determined the numerical values of t(P) up to 0.8 P 0 by assuming 1) that the density of the adsorbate equals that of bulk liquid nitrogen, 2) a cylindrical shape for the pores, the PSD of which was extracted as explained above and 3) that t(P) is independent of the pore radius. For t(P) above 0.8 P 0 , we have chosen to extrapolate the polynomial law which fits the experimental t(P). The isotherm calculated by introducing in Eqs. ( 1) and (2) the t(P) and the PSD shown in the inset of Fig. 1(a) is represented Fig. 3(b).

We will successively discuss predictions (i) and (iii) which are in large discrepancy with the experimental results (Fig. 3(a)); prediction (ii) which is in agreement with experiment and finally, the large shift between the absolute position of the experimental isotherms and the calculated ones.

(i) The adsorption-desorption cycle observed in the case of pores closed at one end is irreversible. The interconnectivity related pore blocking effect described by Mason 9 to explain the origin of the hysteresis phenomenon could not occur because the porous Si layers are composed of non interconnected pores. Moreover, we reject the 'ink bottle' effect, introduced by Everett 10 , as a possible explanation of the isotherm hysteresis loop we observe because of the absence of narrow sections inside the pores.

(iii) The hysteresis loops observed in the porous layer and porous membrane exhibit the same shape of type H2. We have observed, by carrying out adsorption isotherms on porous Si layers with a PSD having mean values varying from 13 nm up to 50 nm, that the position on the pressure axis of the corresponding adsorption and desorption branches are shifted towards P 0 whereas the shape of the hysteresis loops remains of type H2.

Such asymmetry is not expected by the single pore models. Indeed, numerical calculations based on Cohan's model or on DFT 13 predict a symmetrical shape of the hysteresis loop for such large pores in contradiction with what we observe. Furthermore, although porous Si exhibits non-interconnected pores, the hysteresis shape (type H2) we observe is the same as those systematically obtained for highly interconnected mesoporous solids such as controlled pore glass or Vycor 15 . From a theoretical viewpoint, both macroscopic 9 and microscopic 13 approaches described this shape as a signature of the disorder introduced by the interconnectivity between pores. We must conclude that type H2 hysteresis is not characteristic of interconnected porous materials and we thus reject the conclusion made by Ball and Evans according to which a de Boer's shape of hysteresis is a consequence of interconnected network effect.

We believe that the desorption process we observe results from another effect. At the end of the adsorption process, pores are filled by the dense phase and the external surface of the porous layer is covered by a film which connects the pores. This film could play an important role in the desorption process by preventing the concave menisci from passing through each pore at the expected equilibrium pressure, the emptying of the pores being controlled by the tearing of the film at the edges of the pores. The steepness of the desorption branch could be explained if the tearing of the film occurs at a pressure lower than the equilibrium pressures corresponding to the hemispherical menisci of the major proportion of the pores. We note that the desorption of such a more realistic system which consists of a fluid filling and covering pores of finite length with edges at their open ends has never been studied.

Note that the analysis of the desorption branch using BJH method yields a PSD with a width about 0.5 nm which is very far from that of the actual PSD (± 6 nm). It is conventionally admitted 16 that the desorption process in porous materials occurs at coexistence of the high and low dense phases. On such an assumption, a steep desorption branch is interpreted as the signature of the presence of pores of a unique size 17 . Our results obviously show that it is not necessarily the case.

(ii) The sign of the relative positions of the two adsorption branches (P h < P c ) is in agreement with Cohan's model. These experiments reveal two different scenari for the filling of pores open at one or at both ends. Preferential adsorption in the bottom of closed pores is confirmed by analytical calculations performed in an oblique corner formed by the intersection of two planar surfaces 18 . We are currently investigating pore closure effect on adsorption in cylindrical geometry by means of GCMC simulations and DFT calculations.

As shown in Fig. 3(b), the curve which corresponds to the modified Kelvin equation ( 2) is largely outside of the experimental hysteresis loops. For pore sizes as large as in our samples, this equation should account for the equilibrium phase transition. Indeed, microscopic calculations based on DFT performed by Ball and Evans [START_REF] Ball | [END_REF] have shown that, in the case of a distribution of cylindrical pores open at both ends centered on 7.2 nm, which is two times smaller than the mean size of our pores, the coexistence curve is located at a pressure close to that predicted by equation (2). To our knowledge, this is the first time that, for such large pores, such a large discrepancy is found 19 . This may be due to the use of a cylindrical geometry in the calculations which does not account for the actual polygonal shape of pores. GCMC simulations and DFT based calculations are in progress.

In summary, we have studied the adsorption phenomenon in a porous Si medium, a well ordered material with tubular pores without narrow sections and for which we have experimentally proven the absence of interconnection. Type H2 hysteresis are observed whether the pores are open at one or at both ends. These results are in contradiction with Cohan's model and show that the sole observation of type H2 hysteresis is not the signature of the presence of interconnections in porous material as is generally admitted. The steep desorption process cannot be described by the macroscopic Kelvin equation in which we introduce the experimental values of the thickness of adsorbed film prior to the capillary condensation. We note that the application of the widely used BJH method to extract the PSD from this branch is not systematically justified.

Aknowledgments. We thank N. Dupont-Pavlovsky and Professor X. Duval for fruitful discussions and T. McRae for careful reading.

FIG. 1 :

 1 FIG. 1: Bright field TEM images in plane view (a) and crosssection (b) of a porous Si layer with 51% porosity prepared from p + type (∼ 3.10 -3 Ω.cm) [100] Si substrate. Observation axis [100]. Pores (white) are separated by Si walls (black). The pore density is 3.2.10 11 /cm 2 . Inset 1a: PSD. Inset 1b: transmission electronic diffraction pattern.

FIG. 2 :FIG. 3 :

 23 FIG. 2: α spectra from 18 O(p,α) 15 N: Ep = 629 keV. (•) porous Si layer after thermal oxidation in 18 O2 (Ao). (•) part of the layer protected by Al cap (Ac). The arrow indicates the energy of the expected α peak from Ac.


Jean-Yves Fortin
email: jean-yves.fortin@univ-lorraine.fr

Mooyoung Choi
email: mychoi@snu.ac.kr

Distribution of the coalescence times in a system of diffusion-aggregation of particle clusters in one dimension

Keywords: numbers: 05.10.Gg, 05.20.Dd, 05.40.-a, 61.43.Hv

Distribution of the coalescence times in a system of diffusion-aggregation of particle clusters in one dimension

Jean-Yves Fortin, Mooyoung

Introduction

Statistical properties of coalescence and fragmentation processes have been studied in various fields of physics and chemistry [1,2], where, for example, chemical constituents react to form new molecules or fragment into elementary components. Both reactions, aggregation and fragmentation, are governed by stochastic processes. The kinetics of such processes are described by a Smoluchowski equation [3], which is a balance equation for determining the particle distribution in size or mass in the equilibrium state [4,[START_REF] Dammer | Stochastic many-particle systems with irreversible dynamics[END_REF], and can be used for example to determine the size of interstellar dust grains [START_REF] Hayakawa | [END_REF]7].

The dynamical model studied here concerns the properties of aggregating clusters of particles on a one-dimensional chain, where the hard-core condition imposes that only one particle per site is allowed. Only aggregation or coagulation in the absence of fragmentation is studied. These clusters of particles diffuse, collide and merge with other clusters. The merging occurs with probability unity when the distance between two clusters vanishes, according to the rule A∅B → ∅AB or AB∅, where A and B are two clusters of particles merging into a new cluster AB = C with the size given by the sum of the individual sizes of A and B. Also, since the chain is periodic and the process symmetric, there is no internal current. We assume that the merging rate between clusters is independent of the cluster size. Obviously, clusters that are far away from each other will not merge before a characteristic time which depends on their distance and diffusive length. In most of the models described by Smoluchowski equations, such coagulation rates or kernels depend intrinsically on the masses or sizes of the merging components. General descriptions of the physical properties such as mass or size distributions depend indeed on the microscopic details of these kernels [4,8,1,9]. For instance, in the case of a gravitational system, the aggregation rate is proportional to the sum of the masses of the two constituents, eventually with some exponents, whereas in the case of branched polymerization, they are proportional to the product of the molecular weights. Many models can be solved exactly with the help of a generating function [4,10,[START_REF] Connaughton | Kinetics of Cluster-Cluster Aggregation[END_REF], which contains the information on distribution moments. The process of coalescence we consider here is irreversible and produces a single giant cluster within a finite time, which is comparable with gelation phenomena [12,13]. This phenomena has application in the dynamics of random cellular structures such as coarsening of soap froth or bubbles in two dimensions, or evolution of magnetic domains associated with boundaries [14,15]. In this paper, we focus on the universal properties of the coalescence times, in terms of extreme value statistics and scaling behavior in the large-size limit.

The paper is organized as follows: In section 2, we study the local distributions of cluster sizes, the time evolution of which is governed by a stochastic master equation. The generating function for this equation allows us to compute the decreasing concentration of clusters as a function of time, and we compare the results with those obtained via numerical analysis of the stochastic process. In section 3, we consider the distribution of coalescence times in a finite-size system. We show numerically that the limit distribution described by a scaling function of the time divided by the total empty intervals length squared. A diagrammatic time expansion is also proposed in the approximation of independent intervals, and a general formula for the time distribution of an n-interval system is given. A comparison is made with the numerical results of the previous section, as well as modeling with a stochastic system of site diffusing heights in the presence of absorbing conditions.

Master equation for the cluster size distribution

In this section, we introduce the dynamics of clusters defined as a succession of hardcore particles surrounded by at least two empty sites on a periodic chain with N sites of elementary unit length. In particular, we define the probability characterizing the size and position of the cluster on the lattice. Specifically, we write the probability for a cluster of n (≥ 1) particles at time t symbolically in the form

P n (r, t) = Prob(× • r • • • • r+n-1 • •×), ( 1
)
where r is the site location of the first particle of the cluster. Filled and empty circles denote the sites filled with particles and empty sites, respectively while cross symbols (×) represent the site configuration either empty or filled with particles. The system is periodic: P n (r+N, t) = P n (r, t). The number of particles and the number of clusters at time t are given by

N p = n r nP n (r, t) and N c = n r P n (r, t), (2)
respectively. The particle concentration c ≡ N p /N as well as the number N p of particles is conserved. These numbers can be estimated qualitatively via simple arguments. Assuming the simple uncorrelated process for independent particles, we have

P n = N (1 -c) 2 c n and therefore N p = n≥1 nP n = N c. Accordingly, the (mean) number of clusters is equal to N c = n≥1 P n = N c(1 -c),
which thus satisfies a simple semi-circle law, with the maximum cluster number at half filling (c = 1/2).

Stochastic dynamics and the generating function

Let us first consider the general case of fragmentation and coagulation processes.

Clusters fragment or diffuse according to the following rules: At every time step ∆t, a fragment of size m (≤ n) can break out of the cluster from its left boundary with probability l n,m , satisfying n m=1 l n,m = 1. Then it diffuses by one elementary step to the left at rate 1/τ and fuses eventually with another cluster. The same process also occurs on the right side with a fragment of size m and probability r n,m . In order to obtain the Smoluchowski equation, we consider different processes which contribute to the change of probability P n (r, t) after a small time interval ∆t and note that ∂P n /∂t expresses the rate of change of the probability in the limit ∆t → 0. The first process (a) is the decay or destruction of the cluster by fragmentation and diffusion:

∂P n ∂t a = -Prob(× •[ r • • • • •] m • • • • r+n-1 • •×) -Prob(× • r • • • • • [• • • • r+n-1 •] m • ×), (3)
where the notation [• • • • •] m represents a fragment of size m breaking out from a cluster of larger size n with probability l n,m or r n,m . The right-hand side describes the rate-out process due to such decay of the cluster. The first term can thus be explicitly written as

Prob(× •[ r • • • • •] m • • • • r+n-1 • •×) = 1 τ m l n,m P n (r) = 1 τ P n (r) (4)
and similarly the second term. Another rate-out contribution is provided by process (b), in which a fragment close to the original cluster aggregates from the left or right:

∂P n ∂t b = -Prob(×[• • • • •] m • r • • • • r+n-1 • •×) -Prob(× • r • • • • r+n-1 • • [ • • • • •] m ×). (5)
On the other hand, positive (rate-in) contributions come from three different processes. The first one is the aggregation of clusters, labeled process (c):

∂P n ∂t c = Prob(× r-1 [•] 1 • • • • • r+n-1 • •×) + • • • + Prob(×[ r-1 • • • • •] n-1 • r+n-1 • •×) +Prob(× • r • • • • r+n-2 • • [ •] 1 ×) + • • • + Prob(× • r •• [ • • • • •] n-1 ×). ( 6
)
The next process (d) is breaking of a fragment of size n off from a larger cluster on the right or left, provided that at least two consecutive empty sites are present:

∂P n ∂t d = Prob(×[ r-1 • • • • •] n r+n-1 • • ×) + Prob(× • r • [ • • • • •] n ×). (7)
Finally, the cluster can also originate from the fragmentation of a bigger cluster, which we call process (e):

∂P n ∂t e = Prob(ו [ r-m • • • • •] m • • • r+n-1 • •×) + • • • +Prob(× • r • • • • r+n-1 • [• • • • r+n-1+m • ] m •×). (8)
Summation of all the contributions leads the master equation to take the form

τ ∂P n (r) ∂t = -2P n (r) -P n (r) m≥1 P m (r-m-1) -P n (r) m≥1 P m (r+n+1) + n-1 k=1 P n-k (r+k) m≥k P m (r-m+k-1)r m,k + n-1 k=1 P n-k (r) m≥k P m (r+n-k+1)l m,k
which includes terms linear or quadratic in P n . We first notice that for an asymmetric process with the transition probability l n,m = δ m,1 and r n,m = δ m,n , the dynamics produces clusters of various sizes, the distribution of which is critical for some values of the concentration [16,17,18]. We next consider the symmetric case, where only aggregation of clusters occurs with r n,m = l n,m = δ m,n . In this case, the master equation reduces to

τ ∂P n (r) ∂t = -2P n (r) -P n (r) m≥1 P m (r-m-1) -P n (r) m≥1 P m (r+n+1) + n-1 m=1 P n-m (r+m)P m (r-1) + n-1 m=1 P n-m (r)P m (r+n-m+1) +P n (r-1) 1 - m≥1 P m (r+n) + P n (r+1) 1 - m≥1 P m (r-m) , (10)
and the corresponding out-of-equilibrium dynamics leads in a finite time to an absorbing state consisting of a unique cluster.

To study how the cluster concentration decreases with time, it is convenient to introduce the generating function in the Fourier space:

Ĝ(x, k) = r e -ikr n≥1 x n P n (r) = n≥1 x n Pn (k), (11)
where Pn (k) ≡ r e -ikr P n (r) is the Fourier transform, with the momentum taking the values k = 2πm/N (m = 0, • • • , N -1). It is obvious from this definition that

N c = Ĝ(1, 0) and N p = ∂ x Ĝ(x, 0)| x=1 . (12)
The generating function satisfies a master equation which incorporates all the above processes as well as the one-particle dynamics. Specifically, equation ( 10) is transformed into

τ ∂ ∂t Ĝ(x, k) = -2 (1 -cos k) Ĝ(x, k) -(1 + e -ik ) 1 N k ′ e ik ′ Ĝ(x, k ′ ) Ĝ(e -i(k-k ′ ) , k-k ′ ) +(1 + e ik ) 1 N k ′ e -ik ′ Ĝ(xe i(k-k ′ ) , k ′ ) Ĝ(x, k-k ′ ) -Ĝ(1, k-k ′ ) . ( 13
)
The conservation of the number of particles can be checked directly by computing

∂ t N p = ∂ t ∂ x Ĝ(x, 0)| x=1 = 0.
On the other hand, the number of clusters, N c = Ĝ(1, 0), decreases according to

τ ∂ ∂t Ĝ(1, 0) = - 2 N k ′ e ik ′ Ĝ(1, k ′ ) Ĝ(e ik ′ , -k ′ ), (14)
which depends on the contributions of momentum k ( = 0). In the long-time limit, we can consider the solutions for small momenta with x = 1, and assume that Ĝ(e

-i(k-k ′ ) , k-k ′ ) ≃ Ĝ(1, k-k ′ )
. This simplifies equation ( 13) and yields a closed set of equations:

τ ∂ ∂t Ĝ(1, k) ≃ -k 2 Ĝ(1, k) - 2 N k ′ e ik ′ Ĝ(1, k ′ ) Ĝ(1, k -k ′ ), (15)
which contains a diffusive term proportional to k 2 and a nonlinear term corresponding to the coalescence process.

A mean-field analysis of equation ( 15), employing the ansatz Ĝ(1, k) = N δ k,0 g(t), gives simply g(t) = [g(0) -1 + 2t/τ ] -1 , where g(t) = N c /N is the cluster density. A direct numerical evaluation of equation ( 15) also leads to the same behavior. Under this approximation, the cluster density decreases, after a long time, as t -ν with ν = 1. For comparison, we solve equation (10) numerically, discretizing the time with increment ∆t = 0.02, to obtain the exponent ν = 0.82 for N = 200 (see figure 1). Numerical analysis of equation ( 14) also gives the same scaling with an exponent close to unity.

These results are compared with those obtained via a stochastic method using the following rules: Starting from an initial random configuration of particles, we choose a site randomly. If the site is occupied by a particle, the whole cluster to which the particle belongs moves by one step on the right or left with equal probability. This cluster can merge with a neighboring cluster if they are separated by only one empty site. The time evolution of the cluster density obtained from such stochastic dynamics for 10 4 initial random configurations is displayed in the inset of figure 1. For size N = 500, the decay exponent is estimated as ν = 0.96 in the power-law or long-time regime, which is in agreement with the approximate result discussed above. We also notice that in the case of pure particle coagulation, AA → A∅ or AA → ∅A, where particles coalesce on a single site, the exponent is given by ν = 1/2 instead [19]. In this case the particle concentration decreases and is not constant.

Continuum limit analysis

In this section, we examine the continuum limit in which the lattice constant a → 0 and N → ∞ with the length L 0 = N a finite, by introducing the space coordinate x = ra and cluster size l = na. All possible continuum limit should be compatible with equation (2). Distributions are rescaled according to the rule P n (r)/a 2 → P l (x), and the sums are replaced by integrals. In particular, the conserved number of particles becomes a conservation relation for the total length of the clusters

aN p = r a n a(an) P n (r) a 2 → L c = L 0 0 dx Lc 0 dl P l (x)l, (16)
where L c ≡ lim a→0 aN p is the total length. This continuum limit transforms the problem of discrete clusters into a problem of dynamical intervals of size l with probability P l (x) at location x. The dimensionless ratio c = L c /L 0 (≤ 1) is the fraction of the lattice which is occupied by the intervals (or clusters). In accord, the number of clusters is transformed into

N c = r a n a P n (r) a 2 → N c = L 0 0 dx Lc 0 dl P l (x). ( 17
)
For convenience, we also define the characteristic speed v ≡ a/τ and the diffusion constant D ≡ a 2 /τ , and separate equation (10) into two parts, one with all the linear diffusive terms and the other containing the quadratic terms:

1 a 2 ∂P n (r) ∂t = 1 τ P n (r+1) a 2 -2 P n (r) a 2 + P n (r-1) a 2 +v - P n (r) a 2 Np-n m=1 a P m (r-m-1) a 2 - P n (r) a 2 Np-n m=1 a P m (r+n+1) a 2 + n-1 m=1 a P n-m (r+m) a 2 P m (r-1) a 2 + n-1 m=1 a P n-m (r) a 2 P m (r+n-m+1) a 2 - P n (r-1) a 2 Np-n m=1 a P m (r+n) a 2 -a P n (r+1) a 2 Np-n m=1 a P m (r-m) a 2 . (18)
In the continuum limit, the sum of the linear terms contributes in the second order in a since it corresponds to a diffusive process and is therefore proportional to D. In fact, if v (∝ a) is finite, we have D/v ∝ a → 0 and the sum of the linear terms should be negligible compared with the sum of nonlinear terms. Nevertheless, for generality, we keep the diffusive term and obtain

∂P l (x) ∂t = D ∂ 2 P l (x) ∂x 2 + 2v -P l (x) Lc-l 0 dl ′ P l ′ (x-l ′ ) -P l (x) Lc-l 0 dl ′ P l ′ (x+l) + l 0 dl ′ P l-l ′ (x+l ′ )P l ′ (x) . (19)
This equation includes the diffusive linear term, a term corresponding to negative contributions from clusters aggregating with the original cluster of size l on the right or left, and a positive contribution coming from the aggregation of clusters of size l ′ and ll ′ .

Distribution of coalescence times

In this section, we study the distribution of the coalescence times at which all the particle clusters coalesce, starting from a random initial configuration. If we consider the coalescence of n clusters, there is a finite probability that one of the n intervals between two consecutive clusters vanishes at a given time, and this process continues sequentially until only one cluster or interval of size L = (1c)N survives. This coalescence time is uniquely defined for a finite system, and we expect that it scales like L 2 which is given by the diffusion length formula for a random walk process. Figure 2 represents the plot of the coalescence time distribution rescaled with the quantity L 2 ; it shows clearly that all the curves collapse, regardless of the particle concentration c. The universal curve is found to be non-symmetric with an exponential tail whose decay exponent is approximately 4.95, close to the value π 2 /2 as we will demonstrate below. The problem is, in a set of initially empty intervals with different sizes that increase or decrease by diffusion, to compute the first passage time when only one interval has non-zero size. This problem can be mapped onto a q-state Potts model or a voter model, where particle clusters are represented by domain walls or kinks A in the dual lattice between empty intervals made of holes and labeled by letters (a, b, c, • • •) [20,21], as illustrated in figure 3. Just before two domain walls merge, AA → A∅ or AA → ∅A, the remaining empty hole changes its value according to the state of one of the two neighboring intervals with equal probability. This corresponds to the limit q → ∞ in the Potts model since there is no annihilation process involved in the dynamics [22,23]. Indeed, in these theories, the probability of collision of domain walls with annihilation is equal to 1/(q -1) whereas the collisions with coalescence is equal to (q -2)/(q -1), yielding q = ∞ in our model. Moreover, the coalescence process stops when all the empty sites have the same state value.

Let us consider the case with only two intervals of sizes h 1 and h 2 , which corresponds to two clusters of particles under the periodic conditions imposed for the lattice. When one cluster moves, h 1 and h 2 change by one unit in such a way that h 1 → h 1 ± 1 and h 2 → h 2 ∓ 1, subject to the constraint h 1 + h 2 = L; this reduces the problem to that with a single variable h (≡ h 1 ). A representation in terms of height dynamics is shown in figure 4(a), where the height values (lengths of the intervals or distances between clusters) decrease or increase by one unit transferred to the right or left nearest sites which have non-zero height values. The sites reaching a zero value are removed from the dynamics until one site survives with height L. We now consider the probability P t (h|h 0 ) that the first interval has size h at time t, starting with a given initial condition h 0 (0 < h 0 < L) at t = 0. This probability, abbreviated as P t (h) for simplicity, satisfies a differential equation of diffusion

τ ∂P t (h) ∂t = 1 2 [P t (h+1) + P t (h-1) -2P t (h)] (20)
with absorbing conditions P t (0) = P t (L) = 0. This means that some of the flow of P t (h) is absorbed at the extremities h = 0 or L and therefore the distribution is not normalized. The discrepancy accounts for the probability of the first passage. A spectral basis for this distribution in the Fourier space is given by the set of sine functions, sin(kπh/L) with k = 1, • • • , L -1, and a simple analysis of the above differential equation leads to the expression

P t (h|h 0 ) = 2 L L-1 k=1 sin kπh L sin kπh 0 L exp -1 -cos kπ L t . (21)
The probability P t|h 0 of the first passage when h reaches zero or L relates to the time variation of the integrated distribution [START_REF] Redner | First-Passage Processes[END_REF][START_REF] Nagar | [END_REF]. The integrated value is unity until the first passage time, when it begins to decrease. By the help of the summation formula

h sin(kπh/L) = 1 2 [1 -(-1) k ] sin(kπ/L)[1 -cos(kπ/L)] -1
, this probability is finally expressed as a discrete sum (with t in units of τ ) [26,[START_REF] Andrei N Borodin P S A | Handbook of Brownian Motion -Facts and Formulae 2nd ed Probability and Its Applications[END_REF]:

P t|h 0 = -∂ t h P t (h|h 0 ) = 1 L L-1 k=1 [1 -(-1) k ] sin kπ L sin kπh 0 L exp -1 -cos kπ L t . (22)
In case that L ≫ 1 with h 0 /L = α finite, this distribution has the scaling limit tL 2 → t and P t|h 0 /L 2 → P t|α is given by

P t|α = 2π ∞ k=0 (2k + 1) sin [(2k+1)πα] exp - (2k+1) 2 2 π 2 t . ( 23
)
This expression is normalized and invariant under the transformation α → 1α as expected, since we can replace h 1 by h 2 = Lh 1 by symmetry. The asymptotic value is obtained for k = 0 where P t|α decreases exponentially like e -π 2 t/2 . The decay coefficient π 2 /2 is very close to the one obtained from the numerical results in figure 2. Nevertheless, equation ( 23) does not account for the numerical data of figure 2, where many interacting interval configurations are involved. We would thus like to extend this result to a system of n intervals or equivalently, n coalescing clusters of particles (for large n).

The next case is a system composed of three intervals h 1 , h 2 , and

h 3 = L -h 1 -h 2 .
Only two variables are independent and the equation of diffusion reads

∂P t (h 1 , h 2 ) ∂t = 1 6 [P t (h 1 ±1, h 2 ∓1) + P t (h 1 ±1, h 2 ) + P t (h 1 , h 2 ±1) -6P t (h 1 , h 2 )] , (24)
where factor six takes into account the diffusion of each of the three clusters in the two possible directions, left or right. This equation is valid for h 1 and h 2 strictly inside the domain of existence limited by the triangular perimeter which excludes the values of h 1 and h 2 equal to zero or L as well as the case h 1 + h 2 = L. When the system reaches a partial absorbing configuration where one of the three interval has a zero length (namely, h 1 = 0 or h 2 = 0 or h 3 = 0) and only two clusters survive, we should then apply equation (20) for the remaining variable. Accordingly, there is a hierarchy of diffusion equations which have been studied for an open system in the problem of coalescing walkers [28], with statistical independence of successive interval lengths and partial correlations. In their study the coalescence probability of two walkers is p which corresponds here to p = 1. The general version of equation ( 24) in the case of n intervals with the constraint

h 1 + • • • + h n = L can be written as ∂P (n) t ({h i } 1≤i≤n-1 ) ∂t = n-1 i=1 P (n+1) t (h 1 , • • • , h i , ↓, h i+1 , • • • , h n-1 ) + 1 n ∆ 2 h 1 h 1 -∆ 2 h 1 h 2 + ∆ 2 h 2 h 2 + • • • -∆ 2 h n-2 h n-1 + ∆ 2 h n-1 h n-1 P (n-1) t ({h i }), (25) where P (n)
t ({h i }) is the interval distribution for n clusters. We have also defined the discrete Laplacian operators

∆ hh f h ≡ f h+1 + f h-1 -2f h and mixed derivatives 2∆ hh ′ f h,h ′ ≡ ∆ hh f h,h ′ + ∆ h ′ h ′ f h,h ′ -f h+1,h ′ -1 -f h-1,h ′ +1 + 2f h,h ′ .
The down arrow in the term on the first line symbolizes the discrete derivative operator at h = 0 at the corresponding position and the sum represents the multiple collision rates of two clusters coming from configurations with n + 1 intervals [28]. The mixed derivatives account

t=¡ 1 ¡ 2 h 2 0 h 1 (b) h 1 h 3 =0 h 2 h 5 h 4 h 6 (a) Figure 4. (a)
Modelling of the interval dynamics in terms of heights. A site decreases its height value by one unit (here h 4 → h 4 -1) which is transferred on the nearest site, right or left with equal probability, providing it has a non zero value. As example, site 4 can transfer some of its height (in blue color) to sites 2 or 5, as site 3 which has zero height is non active in the dynamics. (b) Schematic time representation diagram for the coalescence of three particles or intervals under implicit periodic boundary conditions. At time τ 2 the first two clusters coalesce and the red interval vanishes. Subsequently, at final time τ 1 = t the remaining two clusters coalesce, reaching the absorbing state in which the blue interval vanishes.

for the particular diffusion of the intervals, which are not locally independent, as the increase of a given interval leads to the decrease of the neighboring interval on the left or right by the same amount. Instead of solving directly equation [START_REF] Nagar | [END_REF], we consider the diagrammatic expansion illustrated in figure 4(b) for three clusters as an example. Given an n-cluster or interval configuration, we first evaluate the probability that one of the n intervals collapses at time τ n-1 (> 0), leading to an (n-1)-interval state. The process is repeated until only one cluster survives at time t = τ n-1 +τ n-2 +• • •+τ 1 . The distribution of times t then gives the probability of coalescence. We also assume that the intervals are independent but satisfy statistically the global constraint h 1 + • • • + h n = L, leaving n -1 independent variables. We finally average over all the interval size configurations to obtain the approximate distribution in the form:

P (n) tot (t) = t 0 P τ n-1 |h 1 dτ n-1 t-τ n-1 0 P τ n-2 |h 2 dτ n-2 • • • t-τ n-1 -•••τ 3 0 P τ 2 |h n-2 P τ 1 |h n-1 dτ 2 {h i } , (26)
where [• • •] {h i } stands for the average over all the interval configurations. Taking the Laplace transform leads to

P(n) tot (s) ≡ L P (n) t , s = n m=2 P(m) s/(m-1) . (27)
The coefficients 1/(m-1) take into account the fact that there is a probability 1/(m-1) to choose an interval among m -1 independent intervals. It relates to the factor 1/n in the diffusion constant in equation [START_REF] Nagar | [END_REF], which rescales the time for each n-cluster configuration. There is a symmetry coefficient n! included in the normalization factor as we have n possible choices in choosing each coalescent pair of clusters in an n-cluster configuration. The Laplace transform P(n) s can be written explicitly as

P(n) s = 1 L L-1 k=1 [1 -(-1) k ] sin (kπ/L) s + 1 -cos (kπ/L) g n (k), (28)
where g n (k) is the interval size average of the term sin(kπh 0 /L) in equation ( 22) for an n-interval configuration and is defined by

g n (k) = {h i ≥1} δ L, h 1 +•••+hn sin (kπh 1 /L) {h i ≥1} δ L, h 1 +•••+hn . ( 29
)
The summation can be performed through the use of the integral representation of the Kronecker delta function. The number of configurations is thus evaluated as

{h i ≥1} δ L,h 1 +•••+hn = {h i ≥1} π -π dθ 2π e iθ(L-h 1 -•••-hn) = π -π dθ 2π e iθL h≥1 e -iθh n = C 1+ǫ dz 2iπ z L-1 (z -1) n = L-1 n-1 , (30)
where the path of the complex integral follows any circle of radius strictly larger than unity. Applying the same technique to the sine function, we obtain

L-1 n-1 g n (k) = n-2 l=0 L-1 n-2-l (-1) l 2 l+1 sin (1-l) kπ 2L + (l+1) π 2 sin (kπ/2L) l+1 + (-1) k+1 sin (n-1) kπ 2L + (n-1) π 2 2 n-1 sin (kπ/2L) n-1 , (31)
which, as L grows large, approaches a finite limit independent of L:

g n (k) = n-2 l=0 (n -1)! (n -2 -l)! (-1) l (kπ) l+1 sin l+1 2 π + (n -1)! (-1) k+1 (kπ) n-1 sin n-1 2 π . (32)
In figure 5, the inverse Laplace transform of equation ( 27) is plotted for n = 8 and n = 12. The agreement with the scaling curve of figure 2 (we have taken c = 0.5 for reference) is indeed manifested for large t.

We also compare the obtained analytical and numerical results with a model based on the stochastic process of equations ( 24) and [START_REF] Nagar | [END_REF]. Specifically, we consider the diffusion of n sites with non-zero random heights, where n corresponds to the number of intervals, see figure 4(a). A site, say i, is randomly chosen (i = 1, . . . , n) and its height value h i is decreased by one unit, while one of the two randomly chosen nearest neighbor sites with non-zero height has its value increased by the same amount. The sites whose height reaches zero value become inactive and are removed, and the number of active sites decreases. The dynamics of local height diffusion stops when the coalescence time is reached, at which all the sites are empty except one, when only a site with the maximum height remains. For n = 10 and n = 20 sites, the resulting distribution of times converges rapidly to the scaling curve of figure 2, with the chosen constraint n i=1 h i = 100. The average number of clusters in a finite system of concentration c is equal to N c(1-c), with the standard deviation given by √ N c(1-c), and we expect that the universal distribution should be correctly dominated by the configurations of a large number of clusters around this mean value. For example, when c = 0.5 and N = 500, this corresponds to n = 125, and for c = 0.1, n = 45. Both cases give the same universal curve of figure 2 which should be achieved for a relatively large number of intervals. We can interpret the coalescence time distribution P

P (n) tot (t) = δ(t -τ 1 -τ 2 -• • • -τ n-1 ) τ 1 ,•••,τ n-1 = ∞ -∞ dλ 2π e iλt n-1 k=1 e -iλτ k τ k ≃ ∞ -∞ dλ 2π e iλt n m=2 k≥1 a (m) k i λ m-1 + ω k = ∞ -∞ dλ 2π e iλt k 2 ,•••,kn≥1 a (2) k 2 • • • a (n) kn (iλ + ω k 2 ) • • • (i λ n-1 + ω kn ) , (33)
where coefficient a k (s + ω k ) -1 [see equation (28)]. The singularities of the integrand are located at λ = iω k l (l-1) in the upper complex plane. This allows us to compute the integral by means of the residue formula, assuming that the singularities are all simple poles:

P (n) tot (t) ≃ n l=2 k 2 ,•••,kn≥1 (l -1)a (2) k 2 • • • a (n) kn j =l ω k j -ω k l l-1 j-1 e -ω k l (l-1)t . (34)
Asymptotically, the dominant contribution of P (n) tot (t) comes from l = 2 and k 2 = 1, which yields an exponential decay:

P (n) tot (t) ≃ c 0 (n)e -(1-cos(π/L))t ≃ c 0 (n)e -(π 2 /2L 2 )t .
The coefficient c 0 can be extracted from equation (34) in the scaling limit t → tL 2 :

c 0 (n) = L 2 a (2) 1 n m=3 L-1 k=1 a (m) k ω k -ω 1 (m -1) -1 , (35) with a (n)
k = L -1 [1 -(-1) k ] sin(kπ/L)g n (k). In the limit of large L, we have lim L→∞ a (m) k ω k -ω 1 (m -1) -1 = 2k π [1 -(-1) k ] k 2 -(m -1) -1 g m (k), (36)
and coefficient c 0 reduces to

c 0 (n) = 4 n m=3 ∞ k=1 2k π [1 -(-1) k ] k 2 -(m -1) -1 g m (k) , (37)
where g m (k) has been evaluated from equation (32). Since g n (k) ≃ (n -1)/(πk) for k large, the sums over k in equation (37) converge as each term behaves asymptotically like k -2 . We have computed c 0 (n) for several values of the interval number n and plot the result in the inset of figure 5 (black circles). These data points fit well with the red line in the inset described by an exponential function c 0c 1 e -c 2 n with c 0 ≈ 24.41, c 1 ≈ 26.16, and c 2 ≈ 0.147. The dashed line in figure 5 represents the asymptotic fit for the curve of concentration c = 0.5 of figure 2 with L = 250 in the region t ≥ 1, and the intersection at the origin gives c 0 ≃ 24.9, which is close to the analytical result in equation (37). It is therefore concluded that the theory describes accurately the asymptotic regime of coalescent clusters.

Conclusion

Stochastic dynamics of coalescent clusters in a finite-size system displays universal behavior which can be accounted for by a model of statistically independent intervals. This is in turn mapped onto a system of diffusing intervals or heights in one dimension, with absorbing conditions. The probability distribution of the interval size configurations follows a hierarchical descending set of diffusive equations between configurations involving successively n, n-1,..., and 1 intervals, and whose approximate solution can be constructed by considering the time averaged product of successive single absorbing state probabilities. We have then used a Laplace transform to express the global distribution of coalescence times as the product of local Laplace distributions, assuming the independence of intervals but taking into account the finite size constraint for a given number of intervals. The results are corroborated quantitatively by a model of positive heights with absorbing conditions on site, where sites reaching zero height value are removed from the dynamics. In this case the universal curve is recovered numerically in the limit of large numbers of initial intervals. The asymptotic behavior of the universal curve, which follows and exponential, is well accounted for by the theory where the numerical factor and decay coefficient of the exponential can be extracted exactly. The method is less accurate near the origin where local correlations between intervals should be taken into consideration. We finally note that this problem is related to other problems involving the survival time distribution or time to complete coalescence of reactants on different geometries, such as coalescence of particles on fully-connected lattices [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF][START_REF] Turban | [END_REF] where the limit distribution has a similar asymptotic exponential form in e -2t , but with a singular behavior in t -7/2 e -π 2 /4t for t small, with t being the rescaled time proportional to the inverse of the number of sites.

+ m≥0 P

 m≥0 m+n (r-1-m)r m+n,n 1 -k≥1 P k (r+n) + m≥0 P m+n (r+1)l m+n,n 1 -k≥1 P k (r-k) + m≥1 P m+n (r)r m+n,m + m≥1 P m+n (r-m)l m+n,m ,

Figure 1 .

 1 Figure 1. Time evolution of the cluster density N c obtained from the discrete time version of equation(10) with time increment ∆t = 0.02, for systems of size N with concentration c = 0.5. The decay exponent, corresponding to the slope of a dashed line, is found to be ν = 0.82 for N = 200. Inset: Time evolution of the cluster density computed from the stochastic process in the system of size N = 500 with the same particle concentration. Averages over 10 4 random initial configurations are taken, with the time step ∆t = 1/N . The decay exponent (given by the slope of the dashed line) is found to be ν = 0.96.

Figure 2 .

 2 Figure 2. Coalescence time distribution for various values of the particle concentration c in the system of size N = 500. Averages over 10 6 initial configurations have been taken for each concentration value. Observed is the scaling collapse with the asymptotic decay exponent found to be 4.95 (as described by the dashed line).

Figure 3 .

 3 Figure 3. Representation of cluster coalescence in terms of moving domain walls (labeled A) in a q-state Potts model, where all the empty site intervals take different state values (a, b, c, • • •). When two clusters merge (AA → A), the hole (labeled c) changes its value according to the state of its nearest neighbors c → a or c → b (in the present case).

Figure 5 .

 5 Figure 5. Coalescence time distribution. The curves with plus symbols are results from the inverse Laplace transform, equation[START_REF] Andrei N Borodin P S A | Handbook of Brownian Motion -Facts and Formulae 2nd ed Probability and Its Applications[END_REF]. The curves with star and cross symbols are obtained from stochastic numerical analysis of diffusive heights located on each site, with n = 10 and n = 20 sites. Averages are performed with 2 × 10 5 initial samples, and the sum of the heights is 100, which corresponds to (1c)N . The curve of figure2with c = 0.5 is given for comparison, and the dashed line is the fit of this curve in the asymptotic region t ≥ 1. Inset: Plot of the coefficient c 0 (filled circles), see equation (35), versus the number of intervals. The red curve is an interpolation with an exponential function (see text), giving an asymptotic estimate c 0 ≈ 24.4.

  tot (t), using the characteristic function of the survival time interval τ k and writing the definition

(m) k and ω k ≡ 1 -

 1 cos(πk/L) are defined from P

Acknowledgments

This work was supported by the National Research Foundation of Korea through the Brain Pool Program (NRF-2018H1D3A2065321) and also through the Basic Science Research Program (Grant No. 2019R1F1A1046285).