Nizar Bouhlel
email: nizar.bouhlel@insa-rennes.fr

Stéphane Méric
email: stephane.meric@insa-rennes.fr

FMCW Radar System for Transponder Identification

In this paper, we deal with a frequency modulated continuous wave (FMCW) radar used for localizing and tracking targets by frequency evaluation of the received radar beat signal. The radar system achieved with a primary radar (reader) and a secondary radar (transponder) is addressed as super high frequency (SHF) radio frequency identification (RFID). Consequently, considering the transponder as an active target, we achieve an identification application thanks to the shift frequency induced by the transponder. Moreover, the impact of the non-linearity behavior of this transponder on the localization performance, is investigated and a solution is proposed for cancelling non-linear effects.

INTRODUCTION

The frequency-modulated continuous wave (FMCW) radars are systems transmitting a continuous and frequency modulated waveform [START_REF] Stove | Linear FMCW radar techniques[END_REF]. These radars transmit and receive simultaneously and without interruption unlike modulated pulsed radars. The FMCW radars are widely used for industrial and academic with many useful applications [START_REF] Ting | FMCW SAR System for Near-distance Imaging Applications-Practical Considerations and Calibrations[END_REF]. The advantages of the FMCW radars are multiples: excellent performance for near-distance applications, excellent range resolution, no blind distance [START_REF] Stove | Linear FMCW radar techniques[END_REF][START_REF] Mcgregor | Switching system for single antenna operation of an S-band FMCW radar[END_REF], easy implementation, low power consumption, miniaturized system design, simplicity and low cost. Particularly, for near-target detection, the FMCW radar is an alternative to the impulse radar [START_REF] Butler | Benefits of Wide Area Intrusion Detection Systems using FMCW Radar[END_REF] as there no limitation for transmitting pulsewidth and for switching time between transmission an reception. An FMCW radar is used to locate targets [START_REF] Al-Alem | High-Resolution On-Chip S -Band Radar System Using Stretch Processing[END_REF][START_REF] Wang | An FMCW Radar Transceiver Chip for Object Positioning and Human Limb Motion Detection[END_REF][START_REF] Ting | FMCW SAR System for Near-distance Imaging Applications-Practical Considerations and Calibrations[END_REF] since it can detect the ranges and velocities of objects. The FMCW radars are also used in other applications such as snow thickness measurement [START_REF] Galin | Validation of Airborne FMCW Radar Measurements of Snow Thickness Over Sea Ice in Antarctica[END_REF], imaging application -SAR imagery [START_REF] Ting | FMCW SAR System for Near-distance Imaging Applications-Practical Considerations and Calibrations[END_REF]-, terrain displacement monitoring [START_REF] Iglesias | Ground-Based Polarimetric SAR Interferometry for the Monitoring of Terrain Displacement Phenomena[END_REF], as well as automotive applications in collision avoidance [START_REF] Boukari | Six-port FMCW collision avoidance radar sensor configurations[END_REF].

Another way to remotely detect objects is the use of the radio frequency identification (RFID) systems. The RFID allows identification and tracking of persons or objects. Over the last decade, RFID systems have received a great attention due to their simple architecture, smaller dimensions, lowcost, low-power and ease of deployment. The RFID systems operating at HF range are used in several applications such as access control, livestock tracking, etc. With the technological advances of RFID (efficient broadband wireless electronics, integration techniques) it is possible to operate with UHF and SHF range allowing a larger operating distance [START_REF] Ghattas | Performance assessment of FMCW radar processing for transponder identification[END_REF][START_REF] Hazry | Study of integration 2.4 GHz and 5.8 GHz in RFID tag[END_REF]. For it is this progress that new application areas open such as tracking of persons or objects, bio-signal recording [START_REF] Omar | Smart patient management, monitoring and tracking system using radio-frequency identification (RFID) technology[END_REF], logistics monitoring [START_REF] Ulku | The Next Generation In PersonnelPeople Tracking: Active RFID Technology Has Allowed for Enhanced Security and Safety[END_REF] and environmental sensing [START_REF] Namjun | A 5.1-W UHF RFID tag chip integrated with sensors for wireless environmental monitoring[END_REF]. Two types of RFID systems are still used nowadays: passive and active. On the one hand, passive UHF transponders can operate up to a distance of 10 meters. On the other hand, active transponders with their own battery allow higher distances for detection and characterization applications. Therefore, besides the standard data transfer function between reader and transponder, distance information could also be evaluated [START_REF] Heidrich | Local positioning with passive UHF RFID transponders[END_REF]. In the case of active configuration, the connected object is considered as a cooperative target. Also, identification and localization of this kind of cooperative target such as RFID backscatter transponder, is well-known for a long time. Thus, we can combine the RFID concept that considers the target as an active device to radar systems which can evaluate the received radar beat signal for localizing and tracking targets.

In this paper, we deal with the system GESTAR (GEnerator System for TArget Radar identification). This system is a radar system for detecting and localizing cooperative targets in complex environments (urban areas, buildings, forests, snow). The system is composed of a small reader (the primary radar) and a shift frequency identification transponder (the secondary radar) responding to the primary radar only when this latter produces an interrogation signal. Indeed, the primary radar transmits a FMCW that is received by the transponder. The FMCW signal is characterized by a RF carrier frequency, a transmitted power and a frequency sweep. The transponder re-transmits this signal by introducing its own shift frequency. The power transmitted by the transponder is captured by the primary radar, signal processed and the target (transponder) is then precisely located by using a constant false alarm rate (CFAR) detection. Finally, the reader displays on a small screen the position (radial distance), the angular position and the shift frequency of the targets. The particularity and the innovation of GESTAR is the absence of blind distance as well as the great autonomy of the transponder. The performances of this system are characterized by several parameters namely: the range resolution, the shift frequency resolution, the angle estimation accuracy and the maximum range detection capacity.

Two kind of transponders are manufactured and used for this study. The first kind of transponder linearly modulates the re-transmitted signal to the primary radar by adding a shift frequency. The other one operates a non-linear modulation. The primary radar is tested either with a linear, or with a nonlinear transponders. Good results are obtained for the linear transponder in the sense that we correctly detect targets. Nevertheless, the beat signal is modified when the non-linear transponder is used and thus the range-shift frequency cartography is altered as well. This alteration is caused by the presence of several harmonics that are added to the spectrum and caused by the aliasing. Therefore, the identification of the target becomes difficult. We propose in this paper to study the effect of nonlinearity behavior of the transponder in the detection and localization. Moreover, we present a method to remove the undesirable harmonics validated by simulations and measurements.

The paper is organized as follows. Section II introduces the principles of RFID localization based on FMCW signal. The FMCW radar performances are presented in section III. Next, the FMCW radar implementation is presented in section IV. The content of section V is related to the non-linearity behaviors of the transponder and the consequences on the detection and localization. Finally, section VI concludes this paper.

RFID LOCALIZATION BASED ON FMCW SIGNAL

Principle of modulated backscattering

Fig. 1 shows the principle of the RFID localization system using FMCW radar system. The primary

s e (t) FM oscillator s r (t) RFID Transponder (i) (f D ) Low-pass filter R s b (t) s g (t)
Figure 1. Setup of FMCW radar based RFID localization system radar system transmits a signal s e (t) assumed to be a linear frequency modulated signal in time interval [0, T ]. The analytic expression of s e (t) is given as follows

s e (t) = A e exp 2π f 0 t + B 2T t 2 (1)
where A e , f 0 , B, and T denote respectively the amplitude, the RF carrier frequency, the sweep bandwidth, and the sweep duration which can be also considered as the modulation duration. In the case of the RFID transponder located at a distance R from the radar that modulates the carrier frequency through a shift frequency f D , the RFID is considered as a secondary radar. This secondary radar transmits the signal s g (t) given as follows

s g (t) = K A e exp 2π f 0 t - R c -f D t + B 2T t - R c 2 (2)
where the parameter K represents the attenuation coefficient with respect to the radar equation (propagation effect, antenna gains and so on) and c is the velocity of the wave. In the same way, the received signal by the radar s r (t) is an attenuated (K parameter) and delayed version of the transmitted signal with a delay of τ = 2R/c. Therefore, the received signal is given as follows

s r (t) = K A e exp 2π f 0 (t -τ ) -f D t - τ 2 + B 2T (t -τ ) 2 (3)
The principle of an homodyne demodulation of the received signal consists in mixing the signal s e (t) with the received signal s r (t). The expression of the beat signal at the mixer output is s b (t) = s e (t)s * r (t) where (*) denotes the conjugate operator. The delay between the transmitted signal and the received signal results in a constant beat frequency called f b which makes possible to measure the distance between the target (transponder) and the radar (reader). The beat signal is given as follows

s b (t) = K A e exp 2π Bτ T -f D t + φ 0 (4)
where

φ 0 = f 0 τ -B τ 2 /(2 T ) -f D R/c
f b = B T τ = 2 B R c T (5)
and the f b value extracted during each modulation duration T allows the determination of R. On the other hand, in case of a frequency shift f D modulation (f D = 0), a FMCW radar can also measure f D by integrating data from several cycles of wave replication over the long time axis (Fig. 2.(b)). Consequently, the beat signal s b (t) could be described both with respect to the time axis t c and t l [START_REF] Ghattas | Performance assessment of FMCW radar processing for transponder identification[END_REF].

s b (t c , t l ) = K A e e 2πφ 0 e 2πf b tc e 2πf D mT (6)
where φ 0 = φ 0 + Bτ m. To establish the range-shift frequency cartography, the beat signal analysis is carried out in the Fourier domain and also two Fourier transforms need to be applied. A more detailed description of specific Fourier transform is adressed hereafter in the section 4.1. Moreover, we apply a specific signal processing to separate the range information (extracted from f b ) from the shift frequency identification (extracted from f D ) as we explain in the section 3.4.

In case of continuous time Fourier analysis over the t c axis, the first Fourier transform denoted

S c b (f, t l ) = T F [s c b (t c , t l )](f )
, is calculated for τ max < t c < T where the c exposant means the continuous case and τ max corresponds to the delay for a maximum distance. The expression of S c b (f, t l ) is given as follows

S c b (f, t l ) = K A e e 2πφ 0 × (T -τ max ) e 2πf D t l × e -π(f -f b )(T +τmax) sinc (π(f -f b )(T -τ max )) . (7)
It is clear that the maximum of S c b (f, t l ) is reached when f = f b . The application of the second Fourier transform over the t l axis for shift frequency processing where 0 < t l < (M -1)T is given as follows

S c b (f, f ) = K A e e 2πφ 0 × (T -τ max ) e -π(f -f b )(T +τmax) sinc (π(f -f b )(T -τ max )) × (M -1)T e -π(f -f D )(M -1)T sinc π(f -f D )(M -1)T (8)
The expression given by ( 8) which is the two-dimensional Fourier transform of the beat signal makes possible to draw up a range-shift frequency cartography where we can identify targets by their distances and their shift frequency. In case of discrete-time signal description, the beat signal is now sampled. The short time t c is described according to the sampling frequency F e , and the long time t l according to the period repetition frequency F s = 1/T . The discrete-time Fourier transform which is connected to the Fourier transform in the continuous case by the Poisson formula, is given as follows

S b (f, f ) = F e F s ∞ n=-∞ ∞ k=-∞ S c b (f -nF e , f -kF s ) = F e F s (T -τ max )(M -1)T e j2πφ 0 × ∞ n=-∞ e -jπ(f -nFe-f b )(T +τmax) sinc[π(f -nF e -f b )(T -τ max )] × ∞ k=-∞ e -jπ(f -kFs-f D )(M -1)T sinc[π(f -kF s -f D )(M -1)T ]. (9)
The details of calculation are given in appendix B.

Angular analysis

In the case of GESTAR system, the primary radar has 2 receive antennas. Thus, it is possible to measure the angular position of targets. A phase comparison monopulse method [START_REF] Skolnik | Radar Handbook[END_REF] is used to estimate the angular position of the targets in the horizontal plane. This method consists in measuring the phase difference φ observed between the two receiving channels:

∆ Σ = -j tan φ 2 (10)
where ∆ and Σ are the difference and sum respectively of the responses of the target on the two receiving antennas. The measured phase φ is related to the horizontal separation d of the two antennas and the position θ of the target as it is shown in Fig. 3.

Determining φ allows us the estimation of the angular position θ of the target.

FMCW RADAR PERFORMANCES

In our case, the radar performances are set up through range resolution and frequency resolution. These resolutions depends upon sweep duration, sweep bandwidth, and pre-processing of the data. In the following subsections, we define the range and the shift frequency resolutions. We also achieve some simulations corresponding to the range-shift frequency cartography and angular position cartography.

Range resolution

The range resolution is defined as the smallest distance between two targets of the same power that can be separated by the processing [START_REF] Chevalier | Principles of Radar and Sonar Signal Processing[END_REF]. For an homodyne FMCW radar, the range resolution can be represented by [START_REF] Piper | Receiver frequency resolution for range resolution in homodyne FMCW radar[END_REF]

∆R = c T 2 B ∆ f b ( 12
)
where ∆f b is the beat frequency resolution which is proportional to the observation time of the signal processed, ie for the beat signal:

∆f b = 1 T -τ max . (13)
Then, the resolution can be written ∆R = c/(2B). It is clear that the time modulation T does not affect the range resolution as long as it remains much greater than τ max .

Shift frequency resolution

The shift frequency f D induced by the RFID transponder is measured by taking into account the beat signal over M sweep durations [START_REF] Stove | Linear FMCW radar techniques[END_REF]. Assuming the displacement of the RFID transponder less than the range resolution, a discrete Fourier transform (DFT) procedure is carried out on axis time t l to estimate the shift frequency f D . Thus, the shift frequency resolution ∆f D is inversely proportional to the processed signal observation time and

∆f D = 1 M T . ( 14
)
By increasing the time modulation T the shift frequency resolution decreases meaning that this resolution is improved.

System parameters

The specific parameters that characterize our radar system are provided in the Table 1. These are the main parameters of radio transmission, modulation, signal processing, transponder and CFAR detection.

To investigate the radar performance and argue the requested maximum visible range (285 m), a power balance study is required. For this, we need the theoretical received power allowing reliable detection probability of the detection P d 0.9

of targets in a Gaussian noise. Indeed, the reflected power of the RFID transponder, denoted P r at the location R of the primary radar antenna is related to the transmitted power P e and is given by the following radar equation [START_REF] Skolnik | Radar Handbook[END_REF]:

P r = P e G e G r λ 2 σ (4π 3 ) R 4 L (15)
where G e and G r stand for the gain of the transmitting and receiving antennas of the FMCW primary radar, λ = c/f 0 the system wave length in free space, R the distance between the reader (primary radar) and the transponder (secondary radar), σ the radar cross section of the transponder considered as a target and L the loss factor. As the acquired signal can be extremely low, reception amplification stages are needed to operate the data received by the primary radar. The amplification chain introduces noise that is added to the received signal. The average power of the thermal noise in reception is defined by:

P n = K b T 0 B n F (16)
where K b is the Boltzman's constant, T 0 the temperature of the receiving chain, F the noise factor and B n the received bandwidth of the radar. In our case, the average noise power P n is given for B n = 20 MHz, T 0 = 290 K and F = 5 dB, and is equal to -96 dBm. In the presence of Gaussian noise, a signal-to-noise ratio (SN R) of 13 dB is required to obtain a probability of detection, denoted P d of 0.9 with false alarm probability P f a = 10 -3 . The SN R = 13 dB is the value used to estimate the radar [START_REF] Skolnik | Introduction to Radar Systems[END_REF]. This value allows to determine the detection capability of the radar under very favorable conditions (free space).

The processing gain of GESTAR system is a function of the number of analysed points N . Moreover, we apply a Nuttall window as a weighing function on the beat signal defined by [START_REF] Iglesias | Ground-Based Polarimetric SAR Interferometry for the Monitoring of Terrain Displacement Phenomena[END_REF]. The processing gain is weighted by the normalized equivalent noise bandwidth (NENBW) [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF] which is equal to 1.9761 for the Nuttall window. Then, the processing gain is P GR = 10 log 10 (1530/NENBW) = 28.89 dB for a range analysis and is P GD = 10 log 10 (1600/NENBW) = 29.08 dB for the shift frequency analysis. Moreover, the global processing P G is given by

P G = P GR + P GD (17)
The received power P r from the RFID transponder below which the detection probability value is too low, is obtained by P r = 13 -P G + P n = -140.8 dBm ( 18)

From [START_REF] Heidrich | Local positioning with passive UHF RFID transponders[END_REF] we can also deduce the corresponding maximum visible range R mvr . Fig. 4 shows the power balance for maximum visible range. Based on [START_REF] Namjun | A 5.1-W UHF RFID tag chip integrated with sensors for wireless environmental monitoring[END_REF] the shift frequency resolution is equal to ∆f D = 7.8 Hz with 1600 points spread over the interval [0, 12.5 kHz] or [12.5 kHz, 25 kHz]. The resolution of the shift frequency analysis is degraded by the Nutall window used in this study. Then, the resolution at -6 dB is given by 2.627/M T = 20.52 Hz where 2.627 corresponds to the normalized width of the main lobe of the Nutall window. For the range resolution, it is given by ∆R = c/2 B = 2.5 m. As mentioned before, this resolution is degraded by the use of Nutall windows and then the resolution at -6 dB is given by 2.627 × ∆R = 6.567 m. Table 2 shows the performance provided by the radar as function of the parameters of Table 1.

Simulations

First of all, as the shift frequency range is given (see Table 1), an average frequency compensation of 20 kHz is performed regardless of the value of the frequency shift introduced by the transponder.

A maximum error of ±1 kHz is thus committed that corresponds of a ±0.2 meter error on the range distance of the transponder. This error which is much lower than the range resolution, is considered negligible.

We perform simulations with the parameters given in Table 1. We establish the range-shift frequency cartography and the angular position cartography. Fig. 5.a shows the detection of 3 targets (T#1, T#2 and T#3) with radial distances [50,150,200] m for T#1, T#2 and T#3 respectively.

Moreover, the shift frequency values of [20, 20.2, 20.6] kHz characterize T#1, T#2 and T#3 respectively. The targets are located at angles θ = [π/30, π/20, 0] which are measured in the azimuth plane. From Fig. 5, it is clear that the targets are correctly detected and the estimated ranges and shift frequencies correspond to the real ranges and frequencies. The detection level decreases when the

FMCW RADAR IMPLEMENTATION

The radar uses a transmitting antenna and two receiving antennas to allow estimation of the angular position of the targets. Since the transmission and reception are simultaneous, we have maintained high isolation to reach good detection performance. In this section, we present the radar system design with the different components. Fig. 6 shows the block diagram of the various modules that composes the radar system. It mainly consists of a power supply module, a signal processing module, transmission and reception modules, and an antenna system. In the following subsections, we present these modules.

Signal processing and display

The signal processing module is implemented in a field-programmable gate array (FPGA) board which is the core of this module. This board accommodates an ALTERA Cyclone V FPGA, a dual 16-bit digitalto-analog converter, a quad 16-bit analog-to-digital converter, a low power double data rate 2 (LPDDR2) memory, and interface circuits. The FPGA is responsible for generating the radar signal transmitted by the radio module, processing the signal received by the two radio reception chains to detect the presence of the transponders, identifying their shift frequency, locating them by estimating their radial distance and their angular position and display the corresponding information on the organic light-emitting diode (OLED) screen of the radar (Fig. 7). To enable real-time radar analysis, most of the processing Figure 7. The designed radar equipped with an OLED display procedures required for these operations is carried out by dedicated hardware functions implemented in the FPGA and controlled by a 32-bit microprocessor (Nios II System Controller). The algorithms based on simple or regular control structures are achieved through hardware implementation such as the algorithms of distance frequency analysis and shift frequency analysis as well as the estimation of the noise level for the detection threshold determination. Algorithms with a complex control structure benefit from a software implementation, such as the algorithm for detecting targets, estimating their parameters and managing the radar display. This structure makes possible to reduce the time required for the detection of the transponders and the estimation of their parameters to a value of less than 200 ms. The implementation of the spectrum analysis of the beat signal is realized via a chirp Ztransform (CZT) which is the Z-transform of a signal along a spiral contour [START_REF] Rabiner | The chirp z-transform algorithm-a lesson in serendipity[END_REF][START_REF] Rabiner | The chirp z-transform algorithm[END_REF]. The advantage of CZT is to evaluate the beat signal frequency only in a zone of interest [F min , F max ]. The user can define the number of analysis point that makes possible to improve the frequency resolution.

We illustrate in Fig. 8 the different steps for achieving a range-shift frequency cartography from the samples of the beat signal given by ( 6). These samples are collected and positioned into a matrix:

Radio transmitter and receiver

This part consists in selecting components allowing to realize the basic functions of the radio chains (amplification, filtering, transposition of frequency, . . . ). Some basic functions are realized thanks to printed structures optimized in simulation with the Microwave Office simulator. Fig. 9 shows a band pass filter which is optimized by simulation.

A radio module consisting of a transmission board and a two-way reception board is made as exhibited in Fig. 10. The two radio boards are assembled in a nickel-plated aluminum shielding mechanism with a lower cover, a separation shield (Fig. 10.b) and a top cover. The assembled radio module (Fig. 10.d) is in the form of a parallelepiped with dimensions 90×55×25.2 mm and weights 122 g. The transmission board (Fig. 10.a) performs frequency synthesis functions and radio frequency transposition, amplification and filtering of the radar signal provided on two channels in quadrature by the signal processing module. A portion of the transmitted signal is picked up and supplied to the

Antenna

As this paper is not dedicated to the antenna studies, we slightly describe the antenna performances.

The study of antennas firstly consists to select the type of radiating element to achieve the specifications about the volume available for the antennas on the artwork of the radar. The choice is thus made on axial-mode helical antennas with circular polarization capable of offering a directivity and a quality of polarization compatible with the specifications for three channels in a reduced volume. The selected compact antenna device is studied and optimized in simulation using the CST Microwave Studio software. Fig. 11.a shows the compact antenna device with three radiating elements. One element is dedicated to the transmittion and the other twos for the reception. This compact antenna with the radome is measured in an anechoic chamber as shown in the Fig. As a conclusion for this part, we summarize the antenna performances we obtain from the measurements, in Table 3.

Measurement results

Indoor and outdoor measurements are achieved for free-space propagation with the radar (primary and secondary) system as outlined above. We expose range, shift frequency and orientation angle values we take over 1500 measurements in indoor context. The Fig. 12.a shows the distance and the estimated shift frequency. The estimated values of distance and shift frequency which are respectively 15 m and 20.6 kHz match the ground truth. The Fig. 12.b exhibits the measured orientation angle where the ground truth value is 1 degree. It is clear that the variance of this parameter is low and that the average is close to the true value. This fluctuation is justified by the presence of noise that affects the estimation of this parameter. Moreover, the monopulse angle estimation is based on data measures from two antennas and does not take into consideration any prior knowledge about the nature of noise. Finally, the Fig. 12 shows the values of the SNR as a function of the number of measurements. The SNR decrease with respect to the number of measurements is small. Nevertheless, this decrease is explained by the depletion of the transponder batteries.

For outdoor configuration, we achieve measurements with the transponder located at 25 m from the radar and identified by a shift frequency of 20.6 kHz. The Fig. 14 corresponds to the range-shift frequency cartography we obtain thanks to the signal processing described in section 4.1. As depicted in Fig. 14, the output values of this range-shift frequency cartography give a distance of 25 m for the transponder and a shift frequency of 20.6 kHz.

NON-LINEAR EFFECTS

This section describes the non-linearity behavior of the transponder and the consequences on the detection and localization. The range measurement with the non-linear transponder is presented and a solution is given to eliminate undesirable harmonic effects. These effects could be considered as target and could be seen as false alarm.

Non-linear modulation effect in time domain

In a linear approach, the transponder is characterized by the frequency shift f D which is added to the transmitted signal from the reader part of the system as described in the part 2.1. The non-linear modulation we use in this study, is exhibited in Fig. 15. The transponder receives the radar signal s e (t -τ ) and then modulates this received signal by its proper signal g(t). The waveform g(t) is considered as the transponder characteristic. For our study, we propose to modelize the non-linear transponder effect with an alternative periodic signal g(t). The period of g(t) is 1/f D . Fig. 15.b shows g(t) as function of time. Since g(t) is periodic, we can derive it into Fourier series

g(t) = +∞ n=-∞ C n exp(j2πnf D t) (19)
where the complex coefficient C n is given by

C n = f D 1/f D 0 g(t) exp(-j2πnf D t) dt = 2 jπn , (n odd) (20)
For τ max < t = t c < T where t c designates the short time variable (see Fig. 8), the signal received by the reader (primary radar) is given, by

s r (t) = g(t).s e (t -τ ). (21)
Finally, for t = t c + t l where t l = mT (see Fig. 8), the expression of the beat signal becomes

s b (t c , t l ) = g(t c + t l ). exp (j2π(f b t c + φ 0 )) . (22)
Details of calculation are given in appendix C. As mentioned in Section 2.1 (equation ( 9)), the Poisson formula is used to connect the twodimensional Fourier transform S b (f, f ) of the sampled beat signal (according to F e and F s ) and that of the two-dimensional Fourier S c b (f, f ). Using the expression [START_REF] Piper | Receiver frequency resolution for range resolution in homodyne FMCW radar[END_REF], S b (f, f ) can be written as follows

S b (f, f ) = F e F s ∞ n=-∞ ∞ k=-∞ S c b (f -nF e , f -kF s ) (23)
Since an anti-aliasing filter is applied in range frequency, the first term nF e is omitted and the expression can be reduced as follows:

S b (f, f ) = F e F s ∞ k=-∞ S c b (f, f -kF s ) = F e F s e j2πφ 0 (M -1)T (T -τ max ) ∞ k=-∞ ∞ n=-∞ C n × e -jπ(f -kFs-nf D )(M -1)T sinc π(f -kF s -nf D )(M -1)T × e -jπ(f -f b -nf D )(T +τmax) sinc (π(f -f b -nf D )(T -τ max )) . (24)
Details of calculation are given in appendix D. Fig. 16 shows the range-shift frequency cartography as the spectrum magnitude of the function (24) where the variable f represents the beat frequency i.e. the target position, and f the shift frequency. The spectrum magnitude corresponds to the configuration where the transponder is located at the distance 50 m meaning that f b = 250 kHz according to [START_REF] Galin | Validation of Airborne FMCW Radar Measurements of Snow Thickness Over Sea Ice in Antarctica[END_REF], and the shift frequency f D = 19.4 kHz. Moreover, the spectrum has several harmonics due to the aliasing effect according to shift frequencies since we use the interval between 12.5 kHz and 25 kHz. In our application, the cut-off frequency of 3 MHz is used for the anti-aliasing filter, that does not allow aliasing in the [0, 3 MHz] band and does not induce effect to the range response.

Nevertheless, considering Fig. 16, many harmonics are detected along the shift frequency that corresponds to the a non-linear transponder behavior. It is obvious that these harmonics are not connected to any transponder. These harmonics correspond to the maximums reached by the functions sinc{π

(f -f b -nf D )(T -τ max )} and sinc{π(f -kF s -nf D )(M -1)T } for frequencies f = f b + nf D and f = kF s + nf D .
These harmonics increase the false alarm rate that is dramatically useless for a detection system as GESTAR.

Table [START_REF] Butler | Benefits of Wide Area Intrusion Detection Systems using FMCW Radar[END_REF] gives the position of these points as a function of (k, n). In the one hand, the real position of the transponder is that of the position of the fundamental characterized by (k, n) = (0, 1) with range frequency f = f b + f D = 269.4 kHz, meaning that the detected range is 54.01 m, and the shift frequency f = f D = 19.4 kHz. The position along the shift frequency axis is correct but the range value is biased by an additional frequency f D . We recall this bias frequency is compensated by removing 20 kHz from the total beat frequency. In the other hand, the harmonics visible in the magnitude of the spectrum and caused by the non-linear transponder behavior, strongly affect the detection and the location of the real targets because they are considered themselves as targets and then gives erroneous interpretations. The solution to this problem is to suppress these harmonics by calculating their positions. In the following subsection, we present a solution to remove the harmonics.

Simulation of removing undesirable harmonics

To remove the aliasing effects that induce false alarms, we propose to modify the modulation time T . Indeed, from the expression of the bidimensional Fourier transform S b (f, f ) of the beat signal given by the equation ( 24), the position of the harmonics depends on the modulation time T = 1/F s . The idea is to use several values of modulation time T and then to multiply two or more range-shift frequency cartographies corresponding to these different values of T . The fundamental frequency corresponding to the true target is retained by the product and the harmonic frequencies are attenuated. It should be noted that the range resolution remains unchanged since it is given by ∆R = c 2B where B is the transmitted signal bandwidth. However, the shift frequency resolution changes since it is inversely proportional to the observation time of the processed signal ∆f D = 1 M T where M is the repetition number of the modulation time. To avoid a big change in the shift frequency resolution, the used modulation times are very close to each other. In this way, one does not affect too much the shift frequency resolution.

The proposed method is applied to various conditions and some simulations are presented here.

The Fig. 17

Indoor and outdoor measurements

Indoor measurement series are performed with the non-linear transponder of the GESTAR project. We expose some curves describing range, shift frequency and orientation angle measurements. Fig. 18.a shows the estimated range and the estimated shift frequency on indoor measurements that are performed 900 times. The estimated values of range and shift frequency which are respectively 30 m and 19.4 kHz correspond to the true values. Fig. 18.b shows the estimated values of the orientation angle (blue color markers) with the real values of the angle (red color markers) with respect to the distance between the reader and the transponder. We have to note that the angle assessments are close to the true value of the angle for distance values between 10 m and 30 m (indoor configuration). Fig. 18.c shows the values of the SNR according to the number of measurements and with respect to the distance between the reader and the transponder. Also, the SNR decrease is seen as an expected phenomena.

Outdoor measurement series are carried out in an open place as shown in Fig. 19. The transponder is placed directly in front of the radar at different ranges and with different offset angles ranging from 0 radian to 0.3 radian. The shift frequency of the transponder is set at f D = 19.4 kHz. Table 5 describes different configurations for the transponder positions (range and angle) and gives the estimated distances, angles and shift frequency. We also measure the received signal power so that to evaluate the received power with respect to the range and the offset angle which gives the location of the transponder with the relation to the reader one. According to the table, we can note that the estimated shift frequency is the same as the real shift frequency which characterizes the transponder. Moreover, the estimated ranges are close to the ground truth distances between the transponder and the reader. The error rate is low for short distances and increases as distance increases. The angle assess remains reliable both for near-zero orientation angles and for small distances up to 52.5 m. In Fig. 20, we show range-shift frequency cartographies provided by the radar system for three modulation times for two transponder positions: R = 12.5 m and R = 85 m. It is clear that we have the same behavior with the simulations exhibited in Fig. 17 i.e. the same where g(.) is defined according [START_REF] Piper | Receiver frequency resolution for range resolution in homodyne FMCW radar[END_REF]. The last integral is evaluated as follows

Figure 2 .

 2 Figure 2. Beat signal over several modulation duration T with respect to shift frequency f D t l = mT with m = 0, 1, . . . , M -1 where M is the total number of modulation duration (see Fig. 2(b)). On the one hand (f D = 0), we define the beat frequency f b by

Figure 3 .

 3 Figure 3. Angular position estimation

Figure 4 .

 4 Figure 4. Power balance for maximum visible range R mvr = 285 m

Figure 5 .

 5 Figure 5. Example of detection of 3 targets with distances [50, 150, 200] m with respective shift frequencies [20, 20.2, 20.6] kHz. The targets are oriented at angles measured with respect to the vertical θ = [π/30, π/20, 0] radian.

Figure 6 .

 6 Figure 6. Block diagram of the various modules that composes the radar system.

Figure 8 .

 8 Figure 8. Instantaneous frequency over M transmitted ramps each line represents the short time dimension and each column gives the long time dimension. The first CZT distance is applied for each column to form a distance profile matrix. The column number k of this new matrix represents the distance profile of the chirp number k. The second CZT is applied to each line for shift frequency processing. The output of this second CZT provides the radar operator with range-shift frequency cartography from the scene.

Figure 9 .Figure 10 .

 910 Figure 9. Optimized bandpass filter by simulation.

Figure 11 .

 11 Figure 11. The realized compact antenna device with three radiating elements (a) and the front radiation patterns measured in an anechoic chamber, as a function of the orientation angle θ (b).

Figure 12 .

 12 Figure 12. Range and shift frequency indoor measurements (a) and orientation angle indoor measurement curves (b).

Figure 13 .

 13 Figure 13. Evolution of the SNR (dB) with respect to the number of measurements (c).

Figure 14 .

 14 Figure 14. Range-shift frequency cartography established for a single target at distance R = 25 m from the radar with a transponder shift frequency f D = 20.6 kHz.

Figure 15 .

 15 Figure 15. Modeling the behavior of the non-linear transponder.

Figure 16 .

 16 Figure 16. Two-dimensional discrete-time Fourier where sampling frequencies are F e = 20 MHz and F s = 12.5 kHz. The shift frequency of the non-linear transponder is f D = 19.4 kHz. The intervals are [0, 2.5 MHz] for the distance and [19 kHz, 21 kHz] for the shift frequency.

Table 4 .

 4 Coordinates of all points of the range-shift frequency cartography as a function of (k, n) values. The frequency f b = 250 kHz corresponds to a target at a distance R = 50 m. (k, n) (0,1) (-3,3) (-31,21) (-59,39) (-62,41) (-90,59) (-93,61) (-121,79) kF s + nf D (kHz) 19.

Figure 17 .

 17 Figure 17. Simulations of range-shift frequency cartography with different modulation times T = [80, 80.1, 80.2]µs given respectively in (a), (b) and (c). The first and second rows are defined for a distance R = 12.5 m and R = 85 m between transponder and reader, respectively. both configurations f D = 19.4 kHz. For R = 12.5 m, each simulation corresponds to a different time modulation T = {80µs, 80.1µs, 80.2µs}. As shown in Fig. 17.a, .b and .c, the harmonics are presents in all the examples and may be confused with targets. As noted before, the position of these harmonics correspond to the maximums reached by S b (f, f ) for f = f b + nf D and f = kF s + nf D . It is clear that the multiplication of two or three images strongly attenuates the harmonics and retains only the fundamental that corresponds to the real target as it is shown in Fig. 17.d. For R = 85 m, Fig. 17.e, .f and .g show the same behavior but with less apparent harmonics due to the more important magnitude of the noise power.

  .a, .b and .c, the harmonics are presents in all the examples and may be confused with targets. As noted before, the position of these harmonics correspond to the maximums reached by S b (f, f ) for f = f b + nf D and f = kF s + nf D . It is clear that the multiplication of two or three images strongly attenuates the harmonics and retains only the fundamental that corresponds to the real target as it is shown in Fig. 17.d. For R = 85 m, Fig. 17.e, .f and .g show the same behavior but with less apparent harmonics due to the more important magnitude of the noise power.

Figure 18 .Figure 19 .

 1819 Figure 18. Range, shift frequency and orientation angle measurements. (a) Estimation of the distance separating the transponder from the radar and the shift frequency. (b) Angle of orientation of the transponder as a function of its position relative to the transponder. The red color corresponds to the true angle values and the blue color to the estimated values. (c) Evolution of the SNR (dB) as a function of the distance between the transponder and the reader.

Figure 20 .- 1 k=0 1 k=0

 2011 Figure 20. Range-shift frequency cartography provided from real measurements for two configurations of distance R = 12.5 m (first row) and R = 85 m (second row) and with different modulation times T = {80, 80.1, 80.2}µs. The figure (d) is the product of range-shift frequency cartographies given by (a), (b) and (c). The figure (h) is the product of range-shift frequency cartographies given by (e), (f) and (g).

  c + t l )e -j2π(f -f b )tc dt c = +∞ n=-∞ C n e j2πnf D t l T τmax e -j2π(f -f b -nf D )tc dt c = (T -τ max ) +∞ n=-∞ C n e j2πnf D t l e -jπ(f -f b -nf D )(T +τmax) × sinc (π(f -f b -nf D )(T -τ max )) (D2)Then, S c b (f, t l ) can be written according to the aboveS c b (f, t l ) = e j2πφ 0 (T -τ max ) +∞ n=-∞ C n e j2πnf D t l × e -jπ(f -f b -nf D )(T +τmax) sinc (π(f -f b -nf D )(T -τ max )) (D3)The second Fourier Transform of S c b (f, t l ), this time, over t l is given byS c b (f, f ) = (M -1)T 0 S b (f, t l )e -j2πf t l dt l = e j2πφ 0 (T -τ max ) +∞ n=-∞ C n e -jπ(f -f b -nf D )(T +τmax) × sinc (π(f -f b -nf D )(T -τ max )) (M -1)T 0 e -j2π(f -nf D )t l dt l = e j2πφ 0 (T -τ max ) +∞ n=-∞ C n e -jπ(f -f b -nf D )(T +τmax) × sinc (π(f -f b -nf D )(T -τ max )) (M -1)T × e -jπ(f -nf D )(M -1)T sinc π(f -nf D )(M -1)T (D4) Finally, we can deduce (24).

Table 1 .

 1 Simulation Setup Parameters

		Carrier frequency fc	5.8 GHz
	Radio	Transmitter power Pe Antenna gains Ge,Gr	25 mW PIRE 12 dBi
		Receiver noise factor F	5 dB
	Antenna	Type Distance between antenna (d)	3D antenna 0.8λ
		Type	Sawtooth
	Modulation	Frequency Sweep B	60 MHz
		Ramp duration T	80 µs
		Total number M	1600
		Sampling Frequency Fe	20 MHz
		Nb. points (range)	1530 points
	Signal Processing	Nb. points (shift frequency)	1600 points
		Total time M T	128 ms
	Transponder	Radar cross-section σ Shift frequency fD	-10 dBm 2 [19 kHz -21 kHz]
		Noise type	Gaussian
	CFAR	SNR detection probability of false alarm P f a	13 dB 10 -3

Table 2 .

 2 Radar performances achieved

	Detection ability
	SNR	≥ 13 dB
	Maximum visible range	285 m
	Time detection	400 ms
	Range estimation
	Measuring range	[0-500 m]
	Resolution	2.5 m
	Resolution at -6 dB	6.87 m
	Shift frequency estimation
	Resolution	7.8 Hz
	Resolution at -6 dB	20.52 Hz
	Angular estimation
	Measuring range	from -15 deg. to 15 deg.

Table 3 .

 3 Main characteristics of the antenna

	Distance between antenna	0.672 λ
	Transmitted antenna
	Transmitter gain	9.4 dBi
	Transmitter beamwidth (H plane)	47 deg.
	Received antenna
	Receiver gain	9.4 dBi / 8.4 dBi
	Receiver beamwidth	63 deg. / 77 deg.
	Transmitter and receiver isolation	-43 dB / -40 dB

Table 5 .

 5 Comparison between the ground truth values and the assessed parameters (angles and ranges). Measurements of shift frequencies and received signal power are added.

	Ground thruth		Measurements
	Angle Range	Angle	Range	Shift frequency Signal to noise
	θ (rad) R (m) estimation (rad) estimation (m) estimation (kHz)	ratio (dB)
		139.6	0.0001	135	19.4	39.17
	0	115.1 83.4	-0.0002 0.0001	117.5 85	19.4 19.4	37.48 42.23
		52.1	0.0000	52.5	19.4	48.53
	0.1053	115.74 52.9	0.0958 0.1564	117.5 52.5	19.4 19.4	35.95 46.42
	0.1946 117.31	0.0075	120	19.4	28.81
	0.3015 120.54	0.0024	125	19.4	16.74

Shift frequency (kHz)

CONCLUSION

In this paper, a FMCW radar system for RFID identification is presented. The range-shift frequency cartography and the range-angular cartography are carried out from the spectral analysis of the beat signal. We define two kind of transponder operations. We define a linear transponder when the shift frequency is linearly induced by the transmitted signal to the reader system. On contrary, the shift frequency of a non-linear transponder is non-linearly added to the transmitted frequency. Also, the assessment of the performance of the systems are presented. A particular attention is given to the nonlinear transponder. The effects from the non-linear behavior of the transponder on the beat frequency spectrum are studied. Furthermore, a signal processing procedure is proposed to mitigate and remove false alarms. The validation of the solutions are presented and validated based on simulated and real FMCW radar data. In perspective, we project the use of other transmitted waveforms as step frequency continuous wave (SFCW) in order to enhance the signal-to-noise ratio.

ACKNOWLEDGMENT

The authors would like to thank the DGA (Direction Générale de l'Armement) for funding this study.

APPENDIX A. BEAT SIGNAL DESCRIPTION

In the case of radar transponder configuration, we consider the transmitted signal which is described with

and the received signal which have been modified by a shift frequency modulation by the transponder device, is

We have to note that the signal magnitude does not take into account the parameters included in the radar equation. After the homodyne demodulation, the beat signal s b (t) = s e (t)s * r (t) is given following

Let t = t c + mT with m = 0, 1, . . . , M -1, the beat signal s b (t) over several cycles of wave replication is given with respect to the time axis t c and mT

Finally, we obtain

which is the beat signal expressed with [START_REF] Iglesias | Ground-Based Polarimetric SAR Interferometry for the Monitoring of Terrain Displacement Phenomena[END_REF].

APPENDIX B. 2D FOURIER TRANSFORM OF THE BEAT SIGNAL

The bi-dimensional Fourier transform of the continuous beat signal s c b (t c , t l ) is given by

The first integral is given as follows

In the same way


Loïc Turban
email: loic.turban@univ-lorraine.fr

Two-species diffusion-annihilation process on the fully-connected lattice: probability distributions and extreme value statistics

Keywords: reaction-diffusion, random walk, fully-connected lattice, extreme value statistics

We study the two-species diffusion-annihilation process, A + B → Ø, on the fully-connected lattice. Probability distributions for the number of particles and the reaction time are obtained for a finite-size system using a master equation approach. Mean values and variances are deduced from generating functions. When the reaction is far from complete, i.e. for a large number of particles of each species, mean-field theory is exact and the fluctuations are Gaussian. In the scaling limit the reaction time displays extreme-value statistics in the vicinity of the absorbing states. A generalized Gumbel distribution is obtained for unequal initial densities, ρ A > ρ B . For equal or almost equal initial densities, ρ A ≃ ρ B , the fluctuations of the reaction time near the absorbing state are governed by a probability density involving derivatives of ϑ 4 , the Jacobi theta function.

Introduction

In the field of non-equilibrium statistical mechanics, reaction-diffusion processes offer the possibility to study the effects of fluctuations on conceptually very simple model systems like the single-species or the two-species annihilation processes [1][2][START_REF] Ben-Avraham | Diffusion and Reactions in Fractals and Disordered Systems[END_REF][START_REF] Schütz | Exactly solvable models for many-body systems far from equilibrium Phase Transitions and Critical Phenomena[END_REF][START_REF] Ódor | [END_REF][START_REF] Henkel | Non-equilibrium phase transitions: absorbing phase transitions[END_REF][START_REF] Ódor | Universality in non-equilibrium lattice systems[END_REF][START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF][START_REF] Täuber | [END_REF].

In a standard mean-field approximation [10], the bimolecular reaction A + B → Ø displays a t -1 asymptotic decay of the particle densities for equal initial values, ρ A (0) = ρ B (0). For unequal densities, ρ A > ρ B , the approach to the absorbing state, ρ B = 0 and ρ A = ρ A (0) -ρ B (0), is exponential. The mean-field approximation assumes that the system remains homogeneous and ignores the effect of spatial correlations in the distribution of reactants, thus giving a lower bound to the actual particle densities [11].

The relevance in low dimensions of initial concentration fluctuations was pointed out by Ovchinnikov and Zeldovich [12] who found a t -3/4 decay in dimension D = 3 for equal initial densities. This result was soon generalized and a t -D/4 decay was proposed for ρ A (0) = ρ B (0) on the basis of numerical simulations, approximate analytical approaches and scaling arguments [13,14] †. The validity of this asymptotic behaviour was later confirmed by establishing rigorous bounds on the particle density [START_REF] Bramson | [END_REF]18] and through a renormalization group study [19,20].

The slowing down of the process is due to the segregation of A and B particles into A-rich and B-rich domains, at the scale of the diffusion length [14,18,21]. The segregation is a consequence of the initial fluctuations of the densities around their mean values. At long time the reaction is efficient only at the interface between the domains and thus slows down. This effect is relevant below the segregation dimension D seg = 4 at which the t -1 homogeneous mean-field decay is recovered. When generalized to q species [22] the problem has a segregation dimension D seg (q) = 4/(q -1) ≥ 2 [23].

When ρ A (0) > ρ B (0) the density of the minority reactant behaves asymptotically as

ρ B (t) ∼ e -λ D G D (t) , (1.1)
with [START_REF] Bramson | [END_REF]18]

G D (t) =      √ t , D = 1 t/ ln t , D = 2 t , D = 3 (1.2)
Note that the upper critical dimension, as for the single-species process, is D c = 2 [20,23]. The t -D/4 behaviour can be actually obtained using mean-field rate equations, provided the inhomogeneity of the system is taken into account. Although there is no qualitative change at D c for equal initial densities, the upper critical dimension signals itself via logarithmic corrections at D c and a stretched exponential decay below D c for unequal initial densities. The two-species annihilation process has potential applications in different domains. It can be used to model particle-antiparticle annihilation in the early universe [13,24], the kinetics of bimolecular chemical reactions [START_REF] Ovchinnikov | Kinetics of Diffusion Controlled Chemical Processes[END_REF][START_REF] Savara | [END_REF] or electron-hole recombination in irradiated semiconductors [27].

The aim of the present work is to study analytically the kinetics of the two-species reaction-diffusion process on the fully-connected lattice with an emphasis on probability distributions. This is a continuation of previous work on the single-species process [28]. Since the lattice with N sites can only be embedded in a (N -1)-dimensional space, taking the thermodynamic limit requires an infinite-dimensional space and one expects mean-field behaviour. Our purpose is to obtain exact results for the particle density and the reaction time in finite-size systems and to study the extreme-value statistics of the reaction time, in the vicinity of the absorbing state, for both equal or unequal initial densities of the reactants. Table 1. s A and s B are the particle numbers, ρ A and ρ B their densities, N is the number of sites and n = N/2. The last line gives the relations for mean values and variances.

s A = s + d s B = s -d s = sA+sB 2 d = sA-sB 2 ρ A = sA N ρ B = sB N x = s n y = d n ρ A = x+y 2 ρ B = x-y 2 x = ρ A + ρ B y = ρ A -ρ B s A -s A = s -s s B -s B = s -s ∆s 2 A = ∆s 2 ∆s 2 B = ∆s 2
The paper is organized as follows. In section 2 we present the model, its mean-field solution when homogeneity is assumed and give a brief description of our results. In section 3 we study the statistics of the number of particles surviving at a given time, first on a finite system and then in the scaling limit. Section 4 is devoted to a similar study of the reaction time, i.e. the time needed to have a given number of particles remaining. This is followed by the conclusion in section 5. Details of the calculations are given in six appendices.

Model, mean field and main results

Model

We consider the two-species reaction-diffusion process, A + B → Ø, on a fully connected lattice with N sites. Let s A and s B be the number of particles of each type with s A ≥ s B and at most one particle per site. In the following we shall use the variables

s = s A + s B 2 , d = s A -s B 2 . ( 2.1)
Thus d = 0 when the initial densities are equal and s = d when the reaction is complete.

A dictionary giving the relations with standard notations is given in table 1.

The system evolves in time through random sequential updates. An update consists of one or two steps. A first site i is selected at random among the N. When this site is occupied by a particle of type A (B) a second site j is randomly selected among the N. If the destination site is occupied by a particle of type B (A), the two particles annihilate and s → s -1 . In all other cases s is unchanged. d is always conserved. At each update the time t is incremented by 1/N so that t = k/N where k is the number of updates. Note that first selecting a site instead of a particle is vital to keep a constant time increment.

The probabilities for the different events are the following:

• s → s ′ = s -1, with probability: 2 s 2 -d 2 N 2 = s A /N i=A × s B /N j=B + s B /N i=B × s A /N j=A (2.2) • s → s ′ = s, with probability: 1 -2 s 2 -d 2 N 2 = 1 -2s/N i=∅ + s A /N i=A × 1 -s B /N j =B + s B /N i=B × 1 -s A /N j =A (2.3)
Note that, contrary to what occurs on finite-dimensional lattices, the initial distribution of the two species does not matter for the fully-connected lattice.

Mean field solution

In the following we always assume that in the initial state, at t = 0, all the sites are occupied, s A + s B = N and we neglect the effect of spatial fluctuations (D > D seg ). The initial value of s is then n = N/2 according to (2.1). In the scaling limit (s.l.), when N and n → ∞, we introduce the scaled variables

x s.l. = s n = ρ A + ρ B , y s.l. = d n = ρ A -ρ B . (2.4)
where ρ A and ρ B are the particle densities. Note that x is the fraction of occupied sites at t and y is a constant giving the asymptotic value of this fraction when t → ∞.

After a small number of updates, ∆k, according to (2.2) the mean value of s is changed by

∆s = -2 s 2 -d 2 N 2 ∆k (2.5)
where on the right the fluctuations of s around its mean value are neglected. In the scaling limit, ∆s → n dx, ∆k → N dt = 2n dt, yielding

dx dt = - s 2 -d 2 n 2 = -(x 2 -y 2 ) , (2.6)
so that dx x 2 -y 2 = 1 2y dx x -y - dx x + y = -dt . (2.7)
Finally the solution satisfying the initial condition, x = 1 when t = 0, is given by:

x = y 1 + y + (1 -y) e -2yt 1 + y -(1 -y) e -2yt . (2.8)
This yields

ρ A = s A N = x + y 2 = y(1 + y) 1 + y -(1 -y) e -2yt , ρ B = s B N = x -y 2 = y(1 -y) e -2yt 1 + y -(1 -y) e -2yt , (2.9)
for the densities of the two species. When y > 0 the approach to the asymptotic values, ρ A = y and ρ B = 0, is exponential. When y → 0 an algebraic decay is obtained:

ρ A = ρ B = 1 2(t + 1)
.

(2.10) The time needed to reach a given value of x is

t = 1 2y ln (1 -y)(x + y) (1 + y)(x -y) , (2.11)
leading to

t = 1 x -1 , (2.12)
when y → 0.

Main results

The evolution of the system is illustrated in figures 1 and 2 giving respectively the probability distribution S N (s, k) of s at different reaction times and the probability distribution T N (s, k) of the reaction time at different values of s.

As expected for an infinite-dimensional system the mean value of s at time t is in agreement with mean-field theory. Asymptotically, it decays as t -1 when s A = s B and approaches its asymptotic value exponentially when s A > s B . The mean value and the variance both scale as n and the fluctuations of s are Gaussian.

For the statistics of the reaction time t three different regimes are observed with the following results in the scaling limit:

• x > y ≥ 0: The reaction is far from complete with s A = O(n), s B = O(n) and s A ≥ s B .
The mean values of the reaction time are the mean-field ones:

t N = 1 x -1 , y = 0 ; t N = 1 2y ln (1 -y)(x + y) (1 + y)(x -y) , y > 0 . (2.13)
The fluctuations are weak, the variance scaling as n -1 :

∆t 2 N = χ y (x) n ; χ 0 (x) = 1 3x 3 - 1 2x + 1 6 , y = 0 ; χ y (x) = 1 2y 2 x x 2 -y 2 - 1 1-y 2 - 1 4y 1+ 1 y 2 ln (1-y)(x+y) (1+y)(x-y) , y > 0 . (2.14)
The probability density is Gaussian:

T(x, θ) = e -θ 2 /[2χy(x)] 2πχ y (x) , θ = n 1/2 (t -t N ) .
(2.15)

• x = y > 0, s = ny + u:
The reaction is close to completion with unequal numbers of particles (s

A -s B = O(n) and u = s B = O(1)
). The reaction time scales logarithmically with n

t N = 1 2y ln 2ny 1 -y 1 + y + γ -H u , (2.16)
where H u = u j=1 1/j is a harmonic number and γ = 0.577215665 . . . is the Euler constant. The variance is independent of n:

∆t 2 N = 1 4y 2 ζ(2) -H (2)
u .

(2.17)

Here H

(2) u is a generalized harmonic number such that H (m) l = l j=1 1/j m and ζ(2) = π 2 /6. The system displays extreme value statistics. The fluctuations are governed by a generalized Gumbel distribution [29,30], indexed by u:

T ′ (u, θ ′ ) = 1 u! exp -(u + 1)(θ ′ + γ -H u ) -e -(θ ′ +γ-Hu) , θ ′ = 2y(t -t N ) . (2.18) • x = y = 0, s ≥ d: The reaction is close to completion with s A = O(1) ≥ s B ≥ 0.
The reaction time grows as n

t N = n(ζ(2) -H (2) s ) , d = 0 ; t N = n 2d (H s+d -H s-d ) , d > 0 , (2.19)
and the variance as

n 2 ∆t 2 N = n 2 (ζ(4) -H (4) s ) , d = 0 ; ∆t 2 N = n 2 4d 2 2ζ(2)-H (2)
s-d -H

(2)

s+d + 1 d (H s-d -H s+d ) , d > 0 , (2.20)
where ζ(4) = π 4 /90. The fluctuations of the reaction time are even stronger and governed now by derivatives of the Jacobi theta function ϑ 4 :

T ′′ (s, θ ′′ ) = e d 2 θ ′′ (s -d)!(s + d)! s m=0 (m 2 + d/dθ ′′ ) ϑ 4 0, e -θ ′′ , θ ′′ = t n .
(2.21)

Number of surviving particles at a given time

In this section we study the probability distribution S N (s, k) giving the probability to have s + d particles of type A and s -d particles of type B remaining after k updates.

As above we assume that the N sites are initially occupied, s A + s B = N.

Master equation

According to (2.2) and (2.3) the master equation governing the evolution of the system takes the following form

S N (s, k) = 1-2 s 2 -d 2 N 2 S N (s, k-1)+2 (s+1) 2 -d 2 N 2 S N (s+1, k-1) , s = d, . . . , n ,(3.1)
with the boundary condition S N (s > n, k) = 0 and the initial condition S N (s, 0) = δ s,n . In (3.1) the first (second) term on the right gives the probability to be in a state with 2s particles (2s + 2 particles) after k -1 updates and to remain in this state (to have two particles annihilating) at the kth update.

Eigenvalue problem

Let us define the column state vector |S N (k) with components S N (s, k), s = d, . . . , n, the master equation (3.1) can be written in matrix form as

|S N (k) = T|S N (k -1)
where the transition matrix T is given by:

T =            1 2 (d+1) 2 -d 2 N 2 0 0 0 0 0 1-2 (d+1) 2 -d 2 N 2 2 (d+2) 2 -d 2 N 2 0 0 0 . . . . . . 0 0 0 1-2 s 2 -d 2 N 2 2 (s+1) 2 -d 2 N 2 0 . . . . . . 0 0 0 0 0 1-2 n 2 -d 2 N 2            . (3.2)
The eigenvalue equation T|v (r) = λ r |v (r) leads to the linear system

1 -2 s 2 -d 2 N 2 -λ r v (r) s + 2 (s+1) 2 -d 2 N 2 v (r) s+1 = 0 , s = d, . . . , n , (3.3) with v (r) n+1 = 0. It is easy to verify that λ r = 1 -2 r 2 -d 2 N 2 , v (r) s = (-1) r-s v (r) r r-s j=1 (s+j) 2 -d 2 r 2 -(s+j-1) 2 when s < r 0 when s > r , (3.4)
solves the eigenvalue problem (3.3). The solution involves the repeated use of the recursion relation

v (r) s = - (s + 1) 2 -d 2 r 2 -s 2 v (r) s+1 , (3.5)
which follows from (3.3). The value of v (r) r , which remains free, will be used to satisfy the initial condition.

Probability distribution S N (s, k)

We look for the initial state vector under the form |S N (0) = n r=d |v (r) which leads to the condition

S N (s, 0) = n r=s v (r) s = δ s,n , s = d, . . . , n , (3.6)
for the components. From the values of v (r) r with r = n, . . . , n -3 (see appendix A) we can infer that the general expression reads:

v (r) r = (n -d)!(n + d)! (n -r)!(n + r)! 2r r -d , r = d, . . . , n . (3.7)
Then, according to (3.4), one obtains:

v (r) s = (-1) r-s 2r(n -d)!(n + d)!(r + s -1)! (n -r)!(n + r)!(s -d)!(s + d)!(r -s)! , s ≤ r . (3.8)
After k updates the state vector |S N (k) is given by

T k |S N (0) = n r=d T k |v (r) = n r=d λ k r |v (r) , (3.9)
which, according to (3.4) and (3.8), gives

S N (s, k) = n r=s λ k r v (r) s = n r=s (-1) r-s 2r(n -d)!(n + d)!(r + s -1)! (n -r)!(n + r)!(s -d)!(s + d)!(r -s)! 1-2 r 2 -d 2 N 2 k (3.10)
for the components.

Mean value and variance when d = 0

Let us define the generating function

S N (w, k) = n s=d sw s S N (s, k) = n r=d r(n -d)!(n + d)! (n -r)!(n + r)! 1-2 r 2 -d 2 N 2 k Ω r,d (w) (3.11)
where

Ω r,d (w) = r s=d (-1) r-s r + s s -d r + d s + d 2sw s r + s .
(3.12) which allows us to evaluate the mean value of s

s N (k) = S N (1, k) = n r=1 2r r j=1 n -j + 1 n + j 1-2 r 2 N 2 k , d = 0 , (3.15)
and its mean-square value:

s 2 N (k) = ∂S N ∂w w=1 = n r=1 2r 3 r j=1 n -j + 1 n + j 1-2 r 2 N 2 k , d = 0 . (3.16)
In the scaling limit (N, n, k → ∞, t = k/N) studied in appendix C, one obtains:

s N (t) n s.l. = 1 t + 1 - 1 6n 1 - 1 (t + 1) 3 , s 2 N (t) n 2 s.l. = 1 (t + 1) 2 - t 2n(t + 1) 4 , d = 0 . (3.17)
In these expressions we kept the sub-leading contributions since the leading ones vanish in the variance given by:

∆s 2 N (t) n s.l. = 1 3(t + 1) - 1 2(t + 1) 3 + 1 6(t + 1) 4 , d = 0 . (3.18)
Thus the fluctuations are small and s is self-averaging. A comparison with finite-size data is shown in figure 3. We were not able to evaluate Ω r,d (w) when d > 0. This case is treated directly in the scaling limit in the next section. Let us assume that in the scaling limit, for any value of d, s N and ∆s 2 N are both growing as n as in (3.17 Furthermore let us write the unknown mean value x as

x = s N n = g y (t) , (3.20)
and define the probability density:

S(σ, t) s.l. = n 1/2 S N (s, k) . (3.21)
Starting from the master equation (3.1) with S N replaced by S, a Taylor expansion of the right-hand-side up to second order in s and k, when re-expressed in terms of the scaled variables, takes the form of an expansion in powers of n -1/2 . The terms independent of n cancel. The terms of order n -1/2 leads to the differential equation

dg y dt = -(g 2 y -y 2 ) , (3.22)
which is the mean-field equation (2.6) so that, according to (2.8),

x = g y (t) = y 1 + y + (1 -y) e -2yt 1 + y -(1 -y) e -2yt , (3.23)
in agreement with (3.17) when y → 0. To the next order, n -1 , one obtains the following partial differential equation:

∂S ∂t = - 1 2 dg y dt + 1 2 dg y dt 2 ∂ 2 S ∂σ 2 + 2g y S + σ ∂S ∂σ .
(3.24) Introducing the reduced variance κ y (t) = ∆s 2 N /n and assuming that S depends on t only through κ y , the partial differential equation (3.24) can be rewritten in the following form:

∂S ∂κ y - 1 2
∂ 2 S ∂σ 2 dκ y dt = - 1 2
dκ y dt +4g y κ y + dg y dt + 1 2 dg y dt 2 ∂ 2 S ∂σ 2 +2g y S+σ ∂S ∂σ +κ y ∂ 2 S ∂σ 2 . ( 3.25)
The left-hand-side and the last bracket on the right-hand-side are both vanishing for the Gaussian density with variance κ y (t), thus we have

S(σ, t) = e -σ 2 /[2κy(t)] 2πκ y (t) , (3.26)
when, according to (3.25), κ y (t) satisfies the first-order differential equation:

dκ y dt + 4g y κ y + dg y dt + 1 2
dg y dt 2 = 0 . (3.27)
The solution is discussed in appendix D and reads: This section is dedicated to the study of T N (s, k), the probability distribution for the number of updates k needed to reach for the first time a total number of surviving particles s A + s B = 2s, as shown in figure 7(a). This probability is related to S N (s + 1, k -1) through:

κ y (t)
T N (s, k) = S N (s + 1, k -1) × 2 (s + 1) 2 -d 2 N 2 . (4.1)
It is given by the product of the probability to be in a state with s ′ = s + 1 after k -1 updates by the probability of the transition s ′ = s + 1 → s ′′ = s at the next update. Making use of (3.10) one obtains:

T N (s, k) = (n -d)!(n + d)! n 2 (s -d)!(s + d)! n r=s+1 (-1) r-s-1 r(r + s)! (n -r)!(n + r)!(r -s -1)! 1-2 r 2 -d 2 N 2 k-1 . (4.2)
The evolution with k is governed by the master equation

T N (s, k) = 1 -2 (s + 1) 2 -d 2 N 2 T N (s, k -1) + 2 (s + 1) 2 -d 2 N 2 T N (s + 1, k -1) , (4.3)
which follows from (3.1) and (4.1).

Generating function

In order to calculate the mean value and the variance of the reaction time t = k/N we introduce the generating function The evolution from s ′ = n to s ′ = s in figure 7(a) proceeds through a succession of steps where the system remains for some time in a state with s ′ A + s ′ B = 2s ′ until two particles annihilate and s ′ → s ′ -1. We associate with such a step the generating function for its lifetime corresponding to the diagrams of figure 7(b):

T N (s, z) = ∞ k=1 z k T N (s, k) . ( 4
L N (s ′ , z) = 1 + z 1 -2 s ′2 -d 2 N 2 + • • • + z l 1 -2 s ′2 -d 2 N 2 l l updates without annihilation + • • • 2z s ′2 -d 2 N 2 annihilation = 2z(s ′2 -d 2 ) N 2 -z [N 2 -2(s ′2 -d 2 )] . (4.5)
The generating function for the reaction time, measured in the number of updates, is obtained as the product:

T N (s, z) = n s ′ =s+1 L N (s ′ , z) = (n -d)!(n + d)! (s -d)!(s + d)! (2z) n-s n s ′ =s+1 {N 2 -z [N 2 -2(s ′2 -d 2 )]} .(4.6)
It is easy to verify on this expression that T N (s, 1) = ∞ k=1 T N (s, k) = 1 so that T N (s, k) is properly normalized.

Mean value and variance

The mean value of the reaction time t = k/N is given by:

t N (s) = 1 N ∞ k=1 kT N (s, k) = 1 N ∂T N ∂z z=1 = n n s ′ =s+1 1 s ′2 -d 2 .
(4.7)

When d = 0 one obtains

t N (s) = n n s ′ =s+1 1 s ′2 = n(H (2) n -H (2) s ) , (4.8)
where H

(m) l = l j=1 1/j m is a generalized harmonic number. When d > 0 one may write

t N (s) = n 2d n s ′ =s+1 1 s ′ -d - 1 s ′ +d = n 2d (H n-d -H s-d -H n+d +H s+d ) , (4.9)
where H l = H

(1) l is a harmonic number. A second derivative gives the mean square value of the reaction time:

t 2 N (s) = 1 N 2 ∞ k=1 k 2 T N (s, k) = 1 N 2 ∂ ∂z z ∂T N ∂z z=1 = n 2 n s ′ =s+1 1 s ′2 -d 2 2 + n 2 n s ′ =s+1 1 (s ′2 -d 2 ) 2 - 1 2 n s ′ =s+1 1 s ′2 -d 2 . (4.10)
Since the first term in the last expression is t N (s)

2 the variance is given by:

∆t 2 N (s) = n 2 n s ′ =v+1 1 (s ′2 -d 2 ) 2 - 1 2 n s ′ =s+1 1 s ′2 -d 2 .
(4.11)

When d = 0 one obtains:

∆t 2 N (s) = n 2 n s ′ =s+1 1 s ′4 - 1 2 n s ′ =s+1 1 s ′2 = n 2 (H (4) n -H (4) s )- 1 2 (H (2) n -H (2) s ) . (4.12)
When d > 0 (4.11) can be rewritten as:

∆t 2 N (s) = n 2 4d 2 n s ′ =s+1 1 (s ′ -d) 2 + 1 (s ′ + d) 2 - 1 4d 1 + n 2 d 2 n s ′ =s+1 1 s ′ -d - 1 s ′ + d = n 2 4d 2 (H (2)
n-d -H

(2)

s-d +H (2)
n+d -H

(2) s+d )-

1 4d 1+ n 2 d 2 (H n-d -H s-d -H n +d +H s+d ) . (4.13)

Scaling limit

In the scaling limit the probability distribution T N (s, k) leads to three different probability densities, depending on the values of x = s/n and y = d/n. We now study these different cases.

4.4.1.

x > y ≥ 0. We first consider the case where both n, s, d → ∞ for fixed values of the ratios x > y ≥ 0. Using in (4.9) the asymptotic expansion for harmonic numbers

H nα = ln n + ln α + γ + O(n -1 ) , (4.14)
where γ is Euler's constant, gives

t N s.l. = h y (x) = 1 2y ln (1 -y)(x + y) (1 + y)(x -y) , (4.15)
which is the mean-field expression (2.11). It reduces to

t N s.l. = 1 x -1 , (4.16)
when y → 0. Using the following expansion for generalized harmonic numbers as well as (4.14), the scaling limit of the variance follows from (4.13) and reads:

H (2) nβ -H (2) nα = nβ j=nα+1 1 j 2 = nβ nα dj j 2 + O(n -2 ) = 1 n 1 α - 1 β + O(n -2 ) (4.17)
∆t 2 N s.l. = χ y (x) n χ y (x) = 1 2y 2 x x 2 -y 2 - 1 1-y 2 - 1 4y 1+ 1 y 2 ln (1-y)(x+y) (1+y)(x-y) . (4.18)
When y → 0 one obtains:

χ 0 (x) = 1 3x 3 - 1 2x + 1 6 . (4.19)
Here too the fluctuations are small and t is a self-averaging variable when x > y.

The dependence on x = ρ A + ρ B of the mean value and the variance of the reaction time t is shown in figure 8 for y = ρ A -ρ B = 0 and figure 9 for y = 1/8. A good collapse of the finite-size data is obtained. The scaling of the variance with n suggests the definition of the following scaleinvariant and centered time variable

θ(s, k) = n 1/2 (t -t N ) , t = k 2n (4.20)
where t N given by (4.15) depends on s through x = s/n. Accordingly we define the probability density as:

T(x, θ) s.l. = 2n 1/2 T N (s, k) . (4.21)
We proceed as in section 3.5 ‡ and solve the master equation ( 4.3) in the scaling limit. ‡ Except that here the expression of t N is already known. We replace T N by T, make a Taylor expansion of T to second order in s and k and rewrite the coefficients and the derivatives in terms of the new scaled variables. In this way an expansion in powers of n -1/2 is obtained. The coefficients of n 0 and n -1/2 vanish identically. The first non-vanishing contribution is coming at order n -1 and leads to the partial differential equation:

∂T ∂x = 1 2(x 2 -y 2 ) 1 2 - 1 x 2 -y 2 ∂ 2 T ∂θ 2 . ( 4.22)
Assuming that T depends on x only through the reduced variance χ y (x) given by (4.18) the partial differential equation transforms into the diffusion equation ∂T ∂χ y = 1 2

∂ 2 T ∂θ 2 . ( 4.23)
Thus the fluctuations are Gaussian:

T(x, θ) = e -θ 2 /[2χy(x)] 2πχ y (x) , x > y ≥ 0 . (4.24)
The Gaussian behaviour of the reaction time is shown for different values of x = ρ A + ρ B in figure 10 for ρ A = ρ B and figure 11 for ρ A -ρ B = 1/8. The finitesize results were obtained by iterating the master equation (3.1), storing the data at each update for a given value of s = nx + 1 and using (4.1).

4.4.2.

x = y > 0, s ≥ d. We study now the case where s = d + u with d = ny, s A = O(n) and s B = u = O(1), i.e. when the reaction is close to completion or complete. The mean reaction time which follows from (4.9)

t N = 1 2y ln 2ny 1 -y 1 + y + γ -H u + O(n -1 ) , (4.25)
has a slow logarithmic growth with n which in the scaling limit yields:

t N ln n = 1 2y + O[(ln n)] -1 ) . (4.26)
The variance in (4.13) gives

∆t 2 N = 1 4y 2 ζ(2) -H (2) u + O ln n n . (4.27)
The fluctuations of the reaction time t are stronger when the reaction is completed or close to completion. The finite-size data collapse for the mean value and the variance as a function of y is shown in figure 12.

In this regime the ratio ∆t 2 N /t N decreases slowly as (ln n) -1 which suggests a new type of statistics. It will be obtained by taking the scaling limit directly on the probability distribution T N (s, k) in (4.2). Since the variance does not depend on n we define a centered time variable as:

θ ′ = 2y(t -t N ) , t = k 2n . (4.28)
The factor 2y is suggested by the form of t N given by (4.25). With this definition one obtains:

k = n y θ ′ + ln 2ny 1 -y 1 + y + γ -H u . (4.29)
Thus one defines the probability density as:

T ′ (u, θ ′ ) s.l. = n y T N (d + u, k) . (4.30)
With the change of summation variable r = j + d + u + 1 (4.2) leads to

n y T N (d + u, k) = 1 u! n-d-u-1 j=0 (-1) j j! A j B j C j D j (1 -E j ) k-1 , (4.31)
where, in the scaling limit with d = ny: Taking into account the expression of k in (4.29) one obtains:

A j = d + j + u + 1 d = 1 + j + u + 1 ny ≃ 1 , (4.32)
B j = (n-d)! (n-d-j -u-1)! = j+u i=0 (n -d -i) ≃ (n-d) j+u+1 , (4.33)
C j = (n+d)! (n+d+j +u+1)! = j+u+1 i=1 (n + d + i) -1 ≃ (n+d) -(j+u+1) , (4.34)
D j = (2d+2u+j +1)! (2d+u)! = j+u+1 i=1 (2d + u + i) ≃ (2d) j+u+1 , (4.35)
E j = (d+j +u+1) 2 -d 2 2n 2 = 2d+j +u+1 2n 2 (j +u+1) ≃ y n (j +u+1) . (4.36)
(1 -E j ) k-1 ≃ 2ny 1 -y 1 + y - (j+u+1)
e -(j+u+1)(θ ′ +γ-Hu) (4.37)

Thus in the scaling limit (4.31) gives

T ′ (u, θ ′ ) = 1 u! e -(u+1)(θ ′ +γ-Hu ∞ j=0 -e -(θ ′ +γ-Hu j j! , (4.38)
which is the generalized Gumbel distribution [29,30]:

T ′ (u, θ ′ ) = 1 u! exp -(u + 1)(θ ′ + γ -H u ) -e -(θ ′ +γ-Hu) . (4.39)
It crosses over to the Gaussian in (4.24) when u ≫ 1. The matching of the two probability densities is studied in appendix E. The collapse of the finite-size data on the generalized Gumbel distribution is shown in figure 13 for different values of u and y = ρ A -ρ B = 1/2. The finite-size data were obtained by iterating (3.1) and using (4.1).

4.4.3.

x = y = 0, s ≥ d. Finally we consider the case when s = O(1) ≥ d. According to (4.8) and (4.12) in the scaling limit the mean value and the variance of the reaction time behave as

t N n s.l. = ζ(2) -H (2) s , ∆t 2 N n 2 s.l. = ζ(4) -H (4) s , (4.40)
when d = 0 whereas (4.9) and (4.13) lead to

t N n s.l. = 1 2d (H s+d -H s-d ) , ∆t 2 N n 2 s.l. = 1 4d 2 2ζ(2)-H (2)
s-d -H

(2) Furthermore, in the scaling limit, one obtains:

s+d + 1 d (H s-d -H s+d ) , (4.41)
1 -2 r 2 -d 2 N 2 k-1 ≃ e -(r 2 -d 2 )θ ′′ , (n -d)!(n + d)! (n -r)!(n + r)! = r-1 j=d n -j n + j + 1 ≃ 1 . (4.45)
Thus the probability density is given by: The collapse of the finite-size data on T ′′ (s, θ ′′ ) is shown at different values of s for s A = s B in figure 16 and for d = (s A -s B )/2 = 4 in figure 17. The finite-size data were obtained by iterating (3.1) and using (4.1). The asymptotic behaviour for θ ′′ ≫ 1 (see figure 18) is governed by the first term in the sum and reads:

T ′′ (s, θ ′′ ) = 2 e d 2 θ ′′ (s -d)!(s + d)! ∞ r=s+1 (-1) r s m=0 (m 2 -r 2 ) e -r 2 θ ′′ . ( 4
T ′′ (s, θ ′′ ) ≃ 2(s + 1)(2s + 1) 2s s + d e -[(s+1) 2 -d 2 ]θ ′′ , θ ′′ ≫ 1 . ( 4

.47)

In order to study the asymptotic behaviour when θ ′′ ≪ 1 it will be convenient to re-express T ′′ (s, θ ′′ ) in terms of Jacobi theta functions. First let us notice that the product in (4.46) vanishes for r ≤ s so that the sum can start at r = 1 instead of s + 1. The product can then be replaced by the differential operator s m=0 (m 2 + d/dθ ′′ ) so that:

∞ r=s+1 (-1) r s m=0 (m 2 -r 2 ) e -r 2 θ ′′ = s m=0 (m 2 + d/dθ ′′ ) ∞ r=1 (-1) r e -r 2 θ ′′ . (4.48) Making use of the identity ( [33] p 463) ∞ r=1 (-1) r q r 2 cos(2rz) = ϑ 4 (z, q) -1 2 , (4.49)
one obtains

T ′′ (s, θ ′′ ) = e d 2 θ ′′ (s -d)!(s + d)! s m=0 (m 2 + d/dθ ′′ ) ϑ 4 0, e -θ ′′ , (4.50)
where the Jacobi theta function can be re-written as ( [START_REF] Whittaker | A Course of Modern Analysis[END_REF] p 475):

ϑ 4 0, e -θ ′′ = π θ ′′ ϑ 2 0, e -π 2 /θ ′′ = 2 π θ ′′ ∞ r=0 e -(r+1/2) 2 π 2 /θ ′′ . ( 4

.51)

Finally the probability density takes the following form:

T ′′ (s, θ ′′ ) = 2e d 2 θ ′′ (s-d)!(s+d)! s m=0 (m 2 +d/dθ ′′ ) π θ ′′ ∞ r=0 e -(r+1/2) 2 π 2 /θ ′′ . (4.52)
The leading contribution when θ ′′ ≪ 1 (see figure 18) comes from the derivatives of the exponential in the first term of the sum so that:

T ′′ (s, θ ′′ ) ≃ 2 √ 2 (s-d)!(s+d)! π 2θ ′′ 2s+5/2 e -π 2 /(4θ ′′ ) , θ ′′ ≪ 1 . (4.53)

Conclusion

The two-species diffusion-annihilation process A + B → Ø has been studied on the fully-connected lattice with size N for either equal or different numbers of particles, s A and s B , in the initial state where s A + s B = N = 2n. Exact probability distributions, S N for s = (s A + s B )/2 at a given time and T N for the reaction time t needed to reach a given number of surviving particles, have been obtained by solving the master equation. The finite-size scaling behaviour of the mean value and the variance of s and t has been determined using a generating function approach.

In the scaling limit, the fluctuations of the number of particles around their mean values, given exactly by mean-field theory, are weak and Gaussian. The statistical properties of the reaction time display three different regimes. When both s A and s B are O(n) the fluctuations of the reaction time, as measured by the ratio R = ∆t 2 N /t N , decay as n -1/2 . The fluctuations are Gaussian and mean-field theory is exact. For unequal initial densities, in the vicinity of the absorbing state when s A = O(n) and s B = O(1), one obtains R = O[(ln n) -1 ] thus the fluctuations are "marginally weak". They are governed by a generalized Gumbel distribution, indexed by s B , which crosses over to the Gaussian density with increasing values of s B . For equal or almost equal initial numbers of particles and when the system is close to the absorbing state, i.e. when s A = O(1) ≥ s B ≥ 0, one obtains R = O(1). The reaction time is then strongly fluctuating. Its probability density involves an alternating infinite series which can be considered as resulting from the applications of a product of s + 1 first-order differential operators to a Jacobi theta function.

A generalized Gumbel distribution has been recently shown to govern the fluctuations of the covering time of the fully-connected lattice, i.e. the time needed for a random walker to visit almost each site at least once [START_REF] Turban | [END_REF] (see also [35]). Although the two problems present some analogies (they share the same value of R and the same probability density), the mean values and the variances scale differently with N, the number of random walkers is constant for the covering time and decreasing for the reaction time.

The extreme value statistics obtained for the reaction time near the absorbing state, for almost equal initial numbers of A and B particles, is quite similar to what was obtained in [28] for the coagulation process A + A → A. Probability densities with nearly the same form of asymptotics govern the behaviour of the squared width of an interface generated by a periodic Brownian motion [36] as well as the maximum height of the 1D Edwards-Wilkinson model [37][38][39][40][41]. i.e. a Gaussian with variance ∆θ 2 = n/(4uy 2 ). This expression has to be compared to (4.24) when x -y = u/n ≪ 1. It is easy to verify that, in this limit, the variance ∆θ 2 = χ y (x) in (4.18) is governed by the first term on the right and reads χ y (x) ≃ n 4uy 2 , (E.8) as expected.

Appendix F. Matching between Jacobian and Gaussian regimes

We were not able to put in evidence this matching at the level of the probability densities. Thus we shall compare the mean values and the variances in the appropriate limits:

x ≪ 1 with x > y ≥ 0 for Gauss and s ≫ 1 with s ≫ d ≥ 0 for Jacobi. For the Gaussian density (2.13)

Figure 1 .

 1 Figure 1. Evolution of the probability distribution of s, S N (s, k), for values of the number of updates k = 2 i , i = 6, . . . , 13 from right to left on a lattice with N = 1024 sites. The reaction is complete when (a) s = d = 0 , (b) s = d = 128. The fluctuations are Gaussian in the scaling limit and maximum for k close to N .

Figure 2 .

 2 Figure 2. Semi-logarithmic plot of the probability distribution T N (s, k) of the number of updates k needed to reach some given values of s = d, d + 2 i , i = 1, . . . , 5 from right to left, for (a) d = 0, (b) d = 64 on a lattice with N = 512 sites. The reaction time isgiven by k/N . The probability distribution evolves from a Gaussian behaviour when s ≫ d to a wide asymmetric extreme-value distribution when the system approaches its absorbing state. Note that the two k scales differ by a factor of 10.

Figure 3 .

 3 Figure 3. Scaling behaviour of (a) the mean value s N and (b) the variance ∆s 2 N of s = (s A + s B )/2 as a function of the time t = k/N . In the initial state ρ A = ρ B = 1/2. The finite-size data for N = 32 (diamond), 64 (square), 128 (circle) were deduced from S N (s, k) given by a numerical iteration of the master equation (3.1). A good collapse on the full lines corresponding to the scaling functions in (3.17) and (3.18) is obtained. The variance is maximum for a number of updates k close to N .

3. 5 .

 5 Scaling limit when d ≥ 0

Figure 4 .

 4 Figure 4. As in figure 3 for y = ρ A -ρ B = 1/8 and N = 64 (diamond), 128 (square), 256 (circle). The scaling functions (full lines) are given in (3.23) and (3.28).

  ) and (3.18) for d = 0. This suggests the introduction, besides the time variable t = k/(2n) and the density y = d/n, of the scaled and centered variable σ(s, k) = s -s N (k)

Figure 5 .

 5 Figure 5. Data collapse for the scaled probability distribution n 1/2 S N (s, k) as a function of σ = n -1/2 (s -s N (k)) at different times t and for increasing lattice sizes, N = 256 (diamond), 512 (square) and 1024 (circle). In the initial state y=0 so that ρ A = ρ B = 1/2. The finite-size data follow from a numerical iteration of the master equation (3.1). The full lines correspond to the Gaussian density (3.26) obtained in the scaling limit. The fluctuations are stronger for t ≃ 1.

  ) 2 e 2yt -(1-y) 2 e -2yt -4y(1+t-y 4 t) 2y(1-y 2 )e -4yt [1+y-(1-y)e -2yt ] 4 . (3.28)In the limit y → 0 (3.18) is recovered. The mean value in(3.23) and the variance are compared to finite-size data in figure4.The scaling behaviour of the probability distribution S N (s, k) at different times is shown in figure5for ρ A = ρ B and figure6for ρ A -ρ B = 1/8.

Figure 6 .

 6 Figure 6. As in figure 5 for y = ρ A -ρ B = 1/8.

4 .

 4 Time required to reach a given number of surviving particles 4.1. Probability distribution

. 4 )Figure 7 .

 47 Figure 7. (a) Evolution of s ′ = (s ′ A + s ′ B )/2 as a function of the number of updates k. Initially all the sites are occupied, s ′ = n. Small circles correspond to the final state of updates where s ′ keeps the same value, bigger circles to a transition s ′ → s ′ -1 when two particles annihilate. (b) Diagrams corresponding to the generating function L N (s ′ , z) in (4.5) for the lifetime of a state with 2s ′ particles.

Figure 8 .

 8 Figure 8. Scaling behaviour of (a) the mean value t N and (b) the variance ∆t 2 N of the time t needed to reach a given value of x = ρ A + ρ B with ρ A = ρ B = 1/2 in the initial state. The finite-size data for N = 64 (diamond), 128 (square), 256 (circle), given by (4.7) and (4.11), collapse on the full lines corresponding to the scaling functions in (4.16) and (4.19).

Figure 9 .

 9 Figure 9. As in figure 8 for ρ A -ρ B = 1/8 and N = 128 (diamond), 256 (square), 512 (circle). The scaling functions are given in (4.15) and (4.18).

Figure 10 .

 10 Figure 10. Data collapse for the scaled probability distribution 2n 1/2 T N (s, k) as a function of θ = n 1/2 (t -t N (s)) at different values of the particle density x = ρ A + ρ B and for increasing lattice sizes, N 64 (diamond), 128 (square) and 256 (circle). In the initial state ρ A = ρ B = 1/2. The full lines correspond to the Gaussian density, T(x, θ) in (4.24), which is obtained in the scaling limit. The fluctuations are growing as x decreases.

Figure 11 .

 11 Figure 11. As in figure 10 for y = ρ A -ρ B = 1/8 with N = 128 (diamond), 256 (square), 512 (circle).

Figure 12 .

 12 Figure 12. Scaling behaviour of (a) the mean value t N and (b) the variance ∆t 2 N of the time t needed to reach a state with y = x = ρ A and u = ρ B = 0 with ρ A = (1+y)/2 and ρ B = (1 -y)/2 in the initial state. The finite-size data for N = 1024 (diamond), 2048 (square), 4096 (circle), 2 16 (triangle), given by (4.7) and (4.11), collapse on the full lines corresponding to the scaling functions in (4.26) and (4.27). The convergence is slow for t N due to a correction to scaling of order (ln n) -1 .

Figure 13 .

 13 Figure 13. Data collapse for the scaled probability distribution nT N (s, k)/y as a function of θ ′ = 2y(t -t N (s)) at different values of u = s B , for y = ρ A -ρ B = 1/2 and increasing lattice sizes, N = 256 (diamond), 512 (square) and 1024 (circle). The full lines correspond to the generalized Gumbel distribution in (4.39), which crosses over to a Gaussian as u increases.

Figure 14 .

 14 Figure 14. Scaling behaviour of (a) the mean value t N and (b) the variance ∆t 2 N of the time t needed to reach a given value of s = (s A + s B )/2 with s A = s B = n in the initial state. The finite-size data for N 64 (diamond), 128 (square), 256 (circle), given by (4.7) and (4.11), given by (4.7) and (4.11), collapse on the full lines corresponding to the scaling functions in (4.40).

Figure 15 .

 15 Figure 15. As in figure 14 with s A = n + 4 and s B = n -4 in the initial state. The scaling functions are now given by (4.41).

Figure 16 .

 16 Figure 16. Data collapse for the scaled probability distribution 2n 2 T N (s, k) as a function of θ ′′ -θ ′′ = (t -t N (s))/n at different values of s = (s A + s B )/2, for s A = s B = n in the initial state and increasing lattice sizes, N = 64 (diamond), 128 (square) and 256 (circle). The full lines correspond to the probability density T ′′ (s, θ ′′ -θ ′′ ) in (4.46), which crosses over to a Gaussian as s increases.

Figure 17 .

 17 Figure 17. As in figure 16 with s A = n + 4 and s B = n -4 in the initial state.

Figure 18 .

 18 Figure 18. Semi-logarithmic plot of T(s, θ ′′ ) in (4.46) for s = 0, 1, 2 with s A = s B = n in the initial state. The dashed lines correspond to the asymptotic behaviour in (4.47) for θ ′′ ≪ 1 and (4.53) for θ ′′ ≫ 1. A similar agreement is obtained when d > 0.

4 )

 4 and the expansion of F a (b) around b 0 gives:F a (b) ≃ -Thus when u ≫ 1 the maximum is amplified, which justifies the approximation, and one obtains:T ′ (u, θ ′ ) ≃ e -uθ ′2 /2 of variables yields T(x, θ) = 2y n 1/2 T ′ (u, θ ′ ) = e -(4uy 2 /n)θ 2 /22π/(4uy 2 /n) , x = y + u n , (E.7)

  and (2.14) yields ≃ ζ(l) -s 1-l /(l -1) for l > 1 andH s±d ≃ ln(s ± d) in (2.19) gives

	t N ≃ for the mean value and 1 x , y = 0 ;		t N ≃	1 2y	ln	x + y x -y	, y > 0 .	(F.1)
	∆t 2 N ≃ for the variance. 1 3nx 3 , y = 0 ; ∆t 2 N ≃	1 4ny 2	1 x -y	+	1 x + y	+	1 y	ln	x-y x+y	, y > 0 , (F.2)
	In the Jacobian regime using the expansions H	(l)
		t N ≃	n s	, d = 0 ;		t N ≃	n 2d	ln	s + d s -d	, d > 0 ,	(F.3) 20)
	leads to
	∆t 2 N ≃	n 2 3s 3 , d = 0 ;	∆t 2 N ≃	n 2 4d 2	1 s -d	+	1 s + d	+	1 d	ln	s -d s + d

s for the mean value, in agreement with (F.1). Finally, given the above expansions, (2., d > 0 , (F.4) for the variance, in agreement with (F.2).

It seems reasonable to conjecture that similar extreme value statistics should govern the reaction time when the system is close to its absorbing state for D > D c = 2 when s A -s B = O(n) in the initial state and for D > D seg = 4 when s A -s B = O(1) in the initial state.

Appendix A. Evaluation of v (r) r According to (3.6), in the initial state:

The recursion relation (3.5) gives

so that:

In the same way

where

, (A.5) so that:

After lengthy but straitforward calculations the same procedure leads to:

. (A.7)

These results suggest the following conjecture

from which (3.7) can be deduced.

Appendix B. Calculation of Ω r,0 (w) and its first derivative at w = 1

According to (3.12) when d = 0

Using the following combinatorial identity for Legendre polynomials [START_REF] Riordan | Combinatorial identities[END_REF] P r (2w

) can be rewritten as:

An integration by parts gives:

When w = 1 one has:

According to the identity [START_REF] Gradshteyn | Tables of Integrals, Series, and Products[END_REF] 1 -1

the integral in (B.5) vanishes and (3.13) is obtained.

The first derivative of (B.4) at w = 1 gives:

The generating function for Legendre polynomials

leads to:

Identifying the coefficients of u r , one obtains Appendix C. Scaling limit of s N /n and s 2 N /n 2 when d = 0

According to (3.15) one may write:

Making use of the following expansions in powers of n -1

one obtains:

(C.5)

In the scaling limit the Euler-Maclaurin summation formula gives

With the change of variable u = r 2 /n (C.5) and (C.6) lead to:

Using ∞ 0 du u a e -bu = a!/b a+1 one obtains s N (k)/n as given in (3.17). The expression of s 2 N (k)/n 2 differs from (C.1) only through a factor u = r 2 /n in the sum over r. Furthermore the correction of order n -1 to the integral in (C.6) now vanishes. Thus we can directly modify the integral in (C.8) to write:

Finally, the integration leads to the expression given in (3.17).

Appendix D. Solution of equation (3.27)

With κ y (t) = α y (t)β y (t) (3.27) leads to

which transforms into a system of two first-order differential equations:

Since g y (t) in (3.23) can be written as

the first equation in (D.2) gives:

For the second equation one has

Where according to (D.1) and (D.3)

The integration in (D.5) is straightforward and gives:

The product of (D.4) and (D.7) leads to: