Student member, IEEE Weike Feng
email: feng.weike.q4@dc.tohoku.ac.jp

Jean-Michel Friedt
email: jean-michel.friedt@femto-st.fr

Grigory Cherniak
email: grigory.chernyak.s6@dc.tohoku.ac.jp

Fellow, IEEE P Motoyuki Sato
email: motoyuki.sato.b3@tohoku.ac.jp

Batch compressive sensing for passive radar range-Doppler map generation

come

1 Abstract-By exploiting the sparsity of the scene containing only a few moving targets, a high-resolution and real-time range-Doppler map generation algorithm for passive bistatic radar is proposed. The proposed algorithm divides the long integration time into multiple short batches, from which a few batches are randomly selected on the basis of compressive sensing theory. A one-dimensional cross-correlation is performed for each selected batch to obtain the range-compressed profile. Mean-value subtraction is then performed to suppress the direct path interference and stationary target reflections. Finally, an extended orthogonal matching pursuit algorithm is proposed for the effective estimation of target Doppler frequency. Practical application of this novel algorithm is examined by the detection of airplanes and ships via two synchronized general-purpose software-defined radio receivers. The results show that the proposed algorithm can achieve an improved resolution and a reduced sidelobe level compared to the conventional algorithms.

Index Terms-Passive bistatic radar, range-Doppler map, software-defined radio, batches algorithm, compressive sensing.

I. INTRODUCTION

ASSIVE bistatic radar (PBR) has attracted considerable attention over the last decades because of its large coverage area, low vulnerability, and reduced electromagnetic pollution to the environment [START_REF] Palmer | An overview of an illuminator of opportunity passive radar research project and its signal processing research directions[END_REF], [START_REF] Howland | Passive bistatic radar systems[END_REF]. Different types of illuminators of opportunity (IO), including frequency modulated broadcast emitters [START_REF] Howland | Fm radio based bistatic radar[END_REF], digital audio broadcasting emitters [START_REF] Coleman | Passive bistatic radar based on target illuminations by digital audio broadcasting[END_REF], spaceborne and Global Navigation Satellite System (GNSS) sources [START_REF] Ma | Maritime moving target indication using passive GNSS-based bistatic radar[END_REF], [START_REF] Li | Passive radar array with lowpower satellite illuminators based on fractional Fourier transform[END_REF], digital communication networks such as Wi-Fi [START_REF] Colone | WiFi-based passive bistatic radar: Data processing schemes and experimental results[END_REF], [START_REF] Feng | Passive radar delay and angle of arrival measurements of multiple acoustic delay lines used as passive sensors[END_REF], and digital terrestrial television broadcasting (DTTB) emitters [START_REF] Gromek | Experimental results of passive SAR imaging using DVB-T illuminators of opportunity[END_REF], [START_REF] Palmer | DVB-T passive radar signal processing[END_REF], have been exploited in PBR to detect and track air and coastal targets. Various signal processing techniques transposed from active radar to passive radar such as space-time adaptive processing (STAP) [START_REF] Dawidowicz | Detection of moving targets with multichannel airborne passive radar[END_REF], [START_REF] Wojaczek | Reciprocal filter-based STAP for passive radar on moving platforms[END_REF], synthetic aperture radar (SAR) imaging [START_REF] Gromek | Experimental results of passive SAR imaging using DVB-T illuminators of opportunity[END_REF], [START_REF] Gromek | Passive SAR imaging using DVB-T illumination for airborne applications[END_REF], and inverse SAR (ISAR) imaging [START_REF] Martorella | Theoretical foundation of passive bistatic ISAR imaging[END_REF], [START_REF] Olivadese | Passive ISAR with DVB-T signals[END_REF] have been developed for specific applications.

In a PBR system, the bistatic range and Doppler frequency of the moving target are estimated by calculating the twodimensional (2D) cross ambiguity function (CAF). Typically, two approaches, namely, the Fourier transform method and the cross-correlation method [START_REF] Howland | Passive bistatic radar systems[END_REF], are used to generate target range-Doppler maps by implementing the fast Fourier transform (FFT). However, for long-term measurements with a high data sampling frequency, the computational burden of these approaches is usually unacceptable, especially when sufficient signal processing gain and fine range-Doppler resolution are desired for a large surveillance area. Therefore, some suboptimal algorithms [START_REF] Stein | Algorithms for ambiguity function processing[END_REF], [START_REF] Langellotti | Comparative study of ambiguity function evaluation algorithms for passive radar[END_REF] have been proposed to obtain real-time processing capability. In particular, the batches algorithm, which is analogous to pulse-Doppler processing in classical active pulse-Doppler radar, has been proposed and validated in [START_REF] Moscardini | Batches algorithm for passive radar: A theoretical analysis[END_REF]. In this algorithm, a long received signal is divided into many short batches. Then, range compression is conducted for each batch by a one-dimensional (1D) crosscorrelation without compensating for the target Doppler shift. Finally, the Fourier transform along each range bin is conducted to estimate the Doppler frequency of the target. The computing time required to generate the target range-Doppler map by the batches algorithm is considerably less than that required by lossless approaches.

When the matched filter theory is used for Doppler estimation, the range-Doppler maps obtained in a limited integration time have a limited resolution and high sidelobe level. The sidelobes produced by the strong targets introduce spurious values in their neighboring Doppler cells, thus leading to wrong position and amplitude estimations of the weak targets. Moreover, FFT will calculate Doppler frequencies up to half of the inverse of the batch duration. When a shorter batch length is used to achieve a lower signal-to-noise ratio (SNR) loss, the maximal calculated Doppler frequency is considerably larger than the Doppler frequency of the target. Therefore, many calculated values of FFT are meaningless and hence ignored. Although a simple way to overcome this is to use the decimation technique before the Fourier transform, this leads to additional computational complexity and SNR loss. Recently, compressive sensing (CS) theory [START_REF] Baraniuk | Compressed sensing[END_REF][START_REF] Ender | On compressed sensing applied to radar[END_REF][START_REF] Weiss | Sparsity-based radar technique[END_REF] has been extensively studied to overcome the sidelobe and resolution problems commonly encountered in the field of radar signal processing, such as ISAR [START_REF] Biao | Sparse coding-inspired highresolution ISAR imaging using multistage compressive sensing[END_REF][START_REF] Qiu | Compressed sensing-based algorithm for passive bistatic ISAR with DVB-T signals[END_REF][START_REF] Kang | Bistatic ISAR imaging and scaling of highly maneuvering target with complex motion via compressive sensing[END_REF], SAR [START_REF] Dong | A novel compressive sensing algorithm for SAR imaging[END_REF][START_REF] Wu | High-resolution passive SAR imaging exploiting structured Bayesian compressive sensing[END_REF][START_REF] Feng | Near range radar imaging based on block sparsity and cross-correlation fusion algorithm[END_REF], STAP [START_REF] Duan | Sparsity-based STAP algorithm with multiple measurement vectors via sparse Bayesian learning strategy for airborne radar[END_REF], [START_REF] Yang | Sparsity-based STAP using alternating direction method with gain/phase errors[END_REF], and passive radar range-Doppler map generation [START_REF] Zheng | Super-resolution delay-Doppler estimation for OFDM passive radar[END_REF], [START_REF] Sevimli | Range-Doppler radar target detection using denoising within the compressive sensing framework[END_REF]. In comparison with conventional approaches, CS theory can achieve higher resolution and better performance. Furthermore, based on the signal sparsity in a specific dictionary and by solving a restricted minimization problem, CS theory can recover the sparse signal from its undersampled measurement with a high probability. Undersampling can reduce the computational complexity and memory usage, making CS theory attractive for practical applications.

In this study, we propose a novel algorithm based on CS to reduce the sidelobe level obtained when using the batches algorithm for CAF calculation in PBR by exploiting the sparsity of the scene. In addition, this algorithm allows processing within a specific Doppler frequency range. For typical applications of PBR, the moving targets are always sparse in the range-Doppler domain after removing direct path interference (DPI) and strong stationary target reflections (such as islands in coastal applications). However, if CS is directly applied to PBR with a long integration time and high sampling frequency, as studied by the previous researchers [START_REF] Weiß | Compressed sensing for passive surveillance radar using dab signals[END_REF], [START_REF] Hadi | Compressed sensing based high-resolution passive bistatic radar[END_REF], the computational complexity may prevent real-time application. Therefore, based on the principle of the batches algorithm, the long integration time is divided into multiple short batches, from which only a few batches are randomly selected and range compression is correspondingly conducted. Finally, we propose an extended orthogonal matching pursuit (EOMP) algorithm to effectively estimate the Doppler frequency of the targets. Compared to the batches algorithm, the proposed algorithm can achieve range-Doppler maps with an improved resolution and a reduced sidelobe level while maintaining the real-time processing capacity. These concepts were experimentally demonstrated using a low-cost PBR system implemented with digital video broadcasting-terrestrial (DVB-T) receivers, as described in detail in the accompanying article describing the hardware setup [START_REF] Feng | Passive bistatic radar using digital video broadcasting-terrestrial receivers as general-purpose software-defined radio receivers[END_REF].

II. SIGNAL MODEL

In a PBR system, the reference signal received by an antenna directly oriented to IOs (e.g., TV towers) can be expressed as

0 ( )= ( ) ( ) ref ref ref ref s t A s t t n t ( 1
)
where Aref is the complex amplitude, s0(t) is the transmitted signal, tref denotes the transmission time from the IOs to the reference antenna, and nref(t) denotes the additive noise.

On the other hand, the received signal from a surveillance antenna, which includes DPI, multipath echoes, and moving target reflections, is given by where χ(τ, f) represents the range-Doppler map of the targets, τ and f are the expected bistatic time delay and Doppler shift of the reflected signal from the target, respectively, (•) denotes the complex conjugate, and Tint denotes the integration time that determines the system Doppler resolution. Notably, to avoid integration losses in (3) [START_REF] Moscardini | Batches algorithm for passive radar: A theoretical analysis[END_REF], the integration time is chosen as Tint = Tref + τmax, where Tref is the length of the reference signal and τmax is the maximal time delay to be processed. In this case, it is assumed that sref(t) = 0 when t (Tref, Tint].

0 0 1 2 0 1 ( ) ( ) ( ) ( ) (
The most straightforward approach to obtain the range-Doppler map of the target is to calculate (3) for each pair of bistatic time delay in (τmin, τmax) and Doppler frequency in (fmin, fmax). However, the computational complexity is normally beyond that manageable in real-time processing, especially for large surveillance areas and fast-moving targets which have large Doppler frequencies. Two simpler methods, namely, the Fourier transform method and the cross-correlation method, can be used. In the first method, the Fourier transform of ssurv(t)•s ref(t -τ) for a given time delay τ is obtained using FFT. In the second method, the cross-correlation between the Doppler-shifted reference signal for a given Doppler frequency f, i.e., sref(t)e -j2πft , and the surveillance signal ssurv(t) is calculated by FFT.

Although these two approaches can reduce the computational cost, two drawbacks limit their applications [START_REF] Moscardini | Batches algorithm for passive radar: A theoretical analysis[END_REF]. First, FFT generates some unused Doppler frequency values in the first method or unused bistatic range values in the second method. Second, if high signal processing gain and range-Doppler resolution are desired, the long integration time and high data sampling frequency may make the calculation of (3) via FFT impractical. To solve these problems, several suboptimal algorithms have been proposed. In particular, for the first problem, based on the prior knowledge of target distance and velocity, the decimation technique is applied before the Fourier transform process. For the second problem, the batches algorithm is proposed, which divides the long integration time into multiple short batches, and this method will be reviewed in the next section.

III. BATCH COMPRESSIVE SENSING ALGORITHM

A. Batches algorithm review

In this subsection, the batches algorithm is briefly reviewed to provide the basis for the proposed algorithm. In the batches algorithm, a dataset collected over a long time is divided into multiple short batches. According to [START_REF] Martorella | Theoretical foundation of passive bistatic ISAR imaging[END_REF], the n-th (n = 1,2,…,N) batch of the reference signal can be expressed as ( , )

( )rect[( / 2)/ ] ref ref b b b S t n s t nT t T T ( 4
)
where Tb is the length of each batch and is assumed to be given by Tb = Tref/N without loss of generality and rect(x) is unity when -1/2 ≤ x ≤ 1/2 and is zero otherwise.

For the surveillance signal, to guarantee the maximum integration gain in each batch, a partially overlapping division approach is employed and the n-th

(n = 1,2,…,N) batch is expressed as max max ( )/2 ( , ) ( )rect b surv surv b b t T S t n s t nT T (5)
Thus, the reference and surveillance signals can be rewritten as

1 1 ( ) ( , ) ( ) ( , ) N ref ref b n N surv surv b n s t S t nT n s t S t nT n (6)
Then, by substituting ( 6) into (3), we get max 2 ( 1)

2 1 0 ( , ) ( , ) ( , ) b b T N j f n T j ft surv ref n f e S tnS t n e d t ( 7
)
Defining the slow time as u = [0,Tb,…,(N -1)Tb] and the nth ambiguity function as

max 2 0 ( , ) ( , ) ( , ) b T j ft n s u r v r e f f S tnS t ne d t (8) (7) can be represented as 2 1 ( , ) ( , ) n N j fu n n f e f (9)
From ( 9), the optimal CAF can be estimated as a weighted sum of the ambiguity functions of all batches. Moreover, if the product of the batch length and the maximum target Doppler frequency is small, then the phase of each sampling in the batch is closer to the phase of the sampling of the central time, i. Based on [START_REF] Wojaczek | Reciprocal filter-based STAP for passive radar on moving platforms[END_REF], the range-Doppler map of the target can be effectively derived using FFT. Specifically, the 1D crosscorrelation can be applied to each batch to obtain the rangecompressed profiles and then the Fourier transform along the slow times can be conducted to estimate the Doppler frequencies of the targets. It should be mentioned that after obtaining the Fourier transform for each range bin, e -jπfTb should be multiplied accordingly. However, if the target Doppler shift can be ignored for each batch, i.e., 12) can be ignored. We have observed that the approximations e -j2πft ≈ e -jπfTb and e -j2πf t ≈ 1 do not show much difference experimentally.

To achieve acceptable approximation precision, the batch length should be short enough to reduce the loss of SNR caused by the batch-wise approximation [START_REF] Martorella | Theoretical foundation of passive bistatic ISAR imaging[END_REF] compared with the optimal process, especially for a fast target with high Doppler frequency. However, the computational complexity of the batches algorithm increases with decreasing batch length (i.e., the number of batches increases). Therefore, a suitable batch length (and hence number of batches) should be considered to ensure a good trade-off between accuracy and computational complexity.

B. Range-Doppler map generation with CS model

Because the matched filter theory is used, conventional algorithms, including the batch algorithm, have a limited range-Doppler resolution owing to the limited frequency bandwidth and integration time. In addition, the Fourier transform-based approach of these algorithms produces high sidelobe levels. The sidelobes of strong targets mask the weaker targets, making processes such as detection and tracking difficult. In general, amplitude weighting methods can be used to reduce the sidelobe level, however, these methods affect the resolution.

CS theory can also be easily applied for PBR [START_REF] Weiß | Compressed sensing for passive surveillance radar using dab signals[END_REF], [START_REF] Hadi | Compressed sensing based high-resolution passive bistatic radar[END_REF]. First, expected time delays in (τmin, τmax) and Doppler frequencies in (fmin, fmax) are discretized into I and J grids. Then, a sparsifying dictionary Ψ can be designed with the ((j -1)I + i)-th column calculated by the reference signal with time delay τi and Doppler frequency fj, i.e., sref(t -τi)e -j2πfjt . Correspondingly, the range-Doppler estimation problem is modeled as surv s n [START_REF] Gromek | Passive SAR imaging using DVB-T illumination for airborne applications[END_REF] where ssurv is the vectorized surveillance signal obtained by stacking ssurv(t) with all samples, χ is the vectorized range-Doppler map, and n is the noise vector. To reduce the complexity and memory used, a random selection approach can be applied, giving un un un surv surv s s n= n [START_REF] Martorella | Theoretical foundation of passive bistatic ISAR imaging[END_REF] where the superscript "un" stands for "undersampling", Φ is the measure matrix used to randomly select rows from Ψ, and the number of rows indicate the selected samplings from ssurv.

Then, based on CS theory and assuming that the range-Doppler vector is sparse, i.e., there are only a few targets in the surveillance area, the solution of ( 14) is given by solving the following optimization problem: Although it has been shown to have non-deterministic polynomial-time hardness, [START_REF] Olivadese | Passive ISAR with DVB-T signals[END_REF] can be approximately solved by many state-of-the-art algorithms, such as OMP algorithm [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF], smoothed l0 norm algorithm [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF], and basis pursuit denoising algorithm [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. However, two problems may make [START_REF] Olivadese | Passive ISAR with DVB-T signals[END_REF] unsuitable for practical application. The first is the considerable computational and memory requirement. Long computing time makes it impossible for real-time application, especially when high sampling rate, long integration time, fine range-Doppler resolution, and large range-Doppler coverage are desired. The second problem, caused by the insufficient sparsity assumption of the range-Doppler map, may occur when DPI is strong and multiple stationary targets are included in the scene. In such a case, the sparsity of the target range-Doppler map is destroyed. Although some algorithms, such as the normalized least mean squares algorithm [START_REF] Ma | A novel two-dimensional sparse-weight nlms filtering scheme for passive bistatic radar[END_REF], [START_REF] Garry | Evaluation of direct signal suppression for passive radar[END_REF], can be employed to suppress DPI and multipath echoes, higher computational complexity is inevitable. In this subsection, we propose a novel CS-based range-Doppler estimation model for PBR to solve these two problems.

First, similar to the process in the batch algorithm, the longtime received signal is divided into multiple short batches. Then, following the Doppler frequency approximation, 1D crosscorrelation is applied to each batch using FFT, which gives us

( ) ( ) ( ) / N tar n n n n N ( 17
)
where the superscript "tar" represents "targets." Therefore, we have a signal matrix consisting of range profiles of all batches, which is defined as

1 m i n m i n 1 m a x m a x ( ) ( ) ( ) ( )
S ( 18
)
where Nbin is the number of range bins. Then, based on ( 12), the batches algorithm can be reformulated to get the target range-Doppler map as

bin Dop N N D C X S
, where ΨD is the Doppler-shifted Fourier transform matrix and is expressed as However, as mentioned above, because of the limited integration time, the Fourier transform-based solution of (20) has a limited resolution and high sidelobe level. For most applications, after DPI and multipath echo suppression, the range-Doppler map is always sparse, i.e., only few elements corresponding to the moving targets in X have significant values. Therefore, if the sparse property of X can be exploited, better resolution and a lower sidelobe level can be obtained. Therefore, based on CS theory, we propose to solve the following restricted optimization problem to estimate the range-Doppler parameters of the targets:

0 T 2 0 0 = min || || , . . || || D F s t X X S X ( 21
)
where

T 0 S S , 0 ( ) D D
, and ||•||F denote the Frobenius norm of a matrix. To further exploit the advantage of CS theory to reduce computational complexity, the batches are randomly selected, resulting in

T 2 0 1 1 = min || || , . . || || F s t X X S X (22)
where 1 1 0

s b i n N N C S S
is the undersampled signal matrix,

0 1 1 s D o p N N D C
denotes the defined sparsifying matrix,

1 s N N C
is the measure matrix used to randomly select batches, and Ns is the number of selected batches. Notably, to save memory usage and reduce computational complexity, undersampling should be performed just after the determination of batches. In such a case, the signal processes mentioned above should be adjusted accordingly.

C. Extended OMP algorithm

In this subsection, we attempt to solve ( 22) effectively and accurately. Although the l0 norm in [START_REF] Biao | Sparse coding-inspired highresolution ISAR imaging using multistage compressive sensing[END_REF] is the best to measure the sparsity of a signal in the absence of noise, its computational burden is very high because of its discontinuity. Therefore, some advanced algorithms should be adopted to solve [START_REF] Biao | Sparse coding-inspired highresolution ISAR imaging using multistage compressive sensing[END_REF]. Compared to [START_REF] Olivadese | Passive ISAR with DVB-T signals[END_REF], the difference in the proposed model is that it is established in the matrix form. For the classical 1D CS problem in [START_REF] Olivadese | Passive ISAR with DVB-T signals[END_REF], three different algorithm types, i.e., greedy algorithms [START_REF] Needell | Cosamp: Iterative signal recovery from incomplete and inaccurate samples[END_REF], [START_REF] Donoho | Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit[END_REF], Bayesian CS algorithms [START_REF] Ji | Bayesian compressed sensing[END_REF], and convex optimization algorithms [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF], have been extensively studied and applied. In general, greedy algorithms are much more computationally efficient than other algorithms and are easier to implement practically. Moreover, among greedy algorithms, the OMP algorithm [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF] is widely used because of its simplicity and efficiency. Therefore, we extend the OMP algorithm, which is normally used to solve sparse vector recovery problems, to the sparse matrix recovery problem in [START_REF] Biao | Sparse coding-inspired highresolution ISAR imaging using multistage compressive sensing[END_REF]. The basic idea of the OMP algorithm is to sequentially find the support set of the measured signal and then project it on the atoms selected from the sparsifying dictionary. Following this basic idea, we propose an extended OMP (EOMP) algorithm, as summarized in Table I.

TABLE I. EOMP ALGORITHM

Input: S1, Ψ1, and iteration number Q. Procedure:

1) Initialization:

1 R S , bin Q N L 0 , 0 Dop bin N N X 0 , 1 q ;
2) Index finding:

H 1 [ , ] argmax | | q q
a b R , where aq and bq are the row index and column index of the biggest element in and(•) H denotes the conjugate transpose; 3) Index set expansion: ( , )

H 1 | | R ,
q q q b a L , and { [ ]} q Ind b K L ,
where Ind{•} denotes the nonzero elements in a vector and Z[u] denotes the u-th column of matrix Z; 4) Estimation:

1 [ ] K , † 1 [ ] q b x S , and 0 [ , ] q b X K x ,
where † ( ) is the pseudo-inverse of a matrix; 5) Residual updating:

1 [ ] [ ] q q b b R S x ; 6) Iteration: 1 q q
, if q Q go back to 2), otherwise, stop.

Output:

T 0

X X .

In each iteration, the atom in Ψ1 that is the most strongly correlated with the residual matrix R (the remaining part of S1) is chosen, and the corresponding row index and column index are recorded in L. Then, for the signal column selected in this iteration, the contribution of all selected atoms is subtracted by projecting this signal column onto the orthogonal subspace spanned by these atoms. The iteration is stopped when the number of iterations equals the maximum iteration number Q. The final estimation of X has non-zero entries corresponding to the index set L, and the values of these entries are obtained by applying the pseudo-inversion process. For the practical implementation for PBR range-Doppler mapping problem, we assume that there are J moving targets with different range-Doppler parameters. In the first iteration, the range and Doppler frequency of the strongest target are estimated by the least squares method. Then, the contribution of this target to the received signal is subtracted. In the next iteration, the secondstrongest target is considered. Finally, all targets are expected to be recovered after Q ≥ J iterations.

D. Signal processing chain and complexity analysis

In summary, the main steps of the proposed range-Doppler map generation algorithm are shown in Fig. 1 and can be described as follows:

1) With the same approach as that shown in [START_REF] Moscardini | Batches algorithm for passive radar: A theoretical analysis[END_REF], divide the reference signal into N consecutive batches and divide the surveillance signal into N partially overlapping batches.

2) Randomly select Ns corresponding batches from the reference and surveillance channels.

3) For each selected batch pair, apply 1D cross-correlation using FFT.

4) Use the mean-value subtraction approach to suppress DPI and stationary target reflections.

5) Given the pre-defined sparsifying dictionary and maximum iteration number, use the proposed EOMP algorithm to generate the range-Doppler map of the targets. Since the real-time processing capacity is considered to be important, the computational complexity of the proposed algorithm is analyzed. The time-consuming parts of the proposed algorithm are the cross-correlation step (item 3 in the algorithm description above) and the EOMP step (item 5). Since FFT can be used for the efficient implementation of the cross-correlation step, the computational burden in terms of complex multiplications for all the selected batches is given by 3NbinNslog2(Nbin). The index finding step dominates the computational consumption of the EOMP algorithm and the complexities of other steps of this algorithm are negligible. For Q iterations, the total complex multiplication number for index finding is given by QNbinNsNDop. Therefore, the complexity of the proposed algorithm can be approximated as [3log2(Nbin)+QNDop]NbinNs. Thus, with an increase in the selected batch number, range bin number, iteration number of EOMP, and Doppler frequency number, the computational load increases. However, thanks to the batch segmentation, undersampling, and the proposed efficient greedy algorithm, real-time range-Doppler generation can be obtained on a general-purpose personal computer (PC), similar to the classical batches algorithm.

IV. EXPERIMENTAL RESULTS

The experimental setup used to collect data is described in [START_REF] Feng | Passive bistatic radar using digital video broadcasting-terrestrial receivers as general-purpose software-defined radio receivers[END_REF]: it mainly consists of two clock-synchronized DVB-T receivers used as general-purpose software-defined radio sources. The sampling rate is 2.048 Msamples/s and data are collected using the GNURadio framework [START_REF] Rondeau | On the GNU Radio ecosystem[END_REF] using two Osmocom sources to display real-time cross-correlation peak values for assessing the quality of the data. As the setup is fully powered by a universal serial bus, it can be easily transported, along with the antennas, for on-field measurements. The experiment setups for landing airplanes and multiple ships are shown in Fig. 2. Fig. 2. Measurements of landing airplanes (left) and multiple ships (right) by the designed PBR system [START_REF] Feng | Passive bistatic radar using digital video broadcasting-terrestrial receivers as general-purpose software-defined radio receivers[END_REF].

A. Fast moving target: landing airplanes

Landing airplanes are selected as the initial target because of their large cross-section and because their large radial velocity v during landing yields Doppler shifts of the order of hundred hertz. Indeed, considering that the reference signal carrier frequency is approximately 500 MHz, the Doppler shift f induced in the carrier frequency fc is given as f = 2fcv/c, where c is the speed of light. When v = 100 m/s, f 330 Hz. This value is the upper bound of the expected Doppler shift introduced by the motion of the airplane in a bistatic configuration assuming that the bistatic angle is null. This value also suggests the integration duration and targeted Doppler coverage.

At the beginning, we show the results obtained by the classical batches algorithm. Fig. 3 Therefore, as shown in Fig. 4 (a), the range-Doppler map can be generated by the Fourier transform in the slow time direction. Furthermore, the range-Doppler map considering 2 14 batches with the same steps of range compression and background subtraction is presented in Fig. 4 (b). By shortening the batch length, a higher SNR can be achieved, which can be demonstrated by the integrated sidelobe level (ISL) of the batch algorithm as a function of batch number, as given in Fig. 5. The ISL in dB scale is defined as

2 ( , ) 2 ( , ) ( ) I , SL 10 log ( ) , i j i j i j i j f f ( 23
)
where Ω is defined as a 3-cell-long and 3-cell-wide window around the target. For short batches, the Doppler shift along the batch length can be ignored. Therefore, as observed in Fig. 5, the sidelobe level is lowered with increasing number of batches, which is consistent with the theoretical analysis. However, the increase in the number of batches induces an increase in the computation time. Furthermore, as shown in Fig. 5, the SNR level is almost constant when the number of batches is greater than 2 14 . Therefore, for the proposed algorithm, 2 14 is selected as the base value. Fig. 6 shows the range-Doppler map obtained by the proposed algorithm with 2 11 batches randomly selected from the 2 14 batches. The iteration number of the proposed EOMP algorithm is set to 10, which is sufficient to estimate the target. Fig. 6 shows that the proposed algorithm can provide a clearer range-Doppler map with only one outlier (caused by the noise), as indicated by the white square. To show the stability of the proposed algorithm with randomly selected undersampled data, 100 random selections of the batches are conducted and the average value is presented. It is experimentally proved that, with 100 trials, the performance (e.g., the ISL) of the proposed algorithm is statistically converged, as shown in Fig. 7 (a).

Taking 100 trials will induce many more outliers than the single trial case, as can be seen from Fig. 7 (b). However, for 100 trials, the target can always be recovered, which demonstrates that no matter how random the selected batches are, the proposed algorithm is robust enough to work well. We note that, in practice, only one-time random selection will be conducted to achieve the real-time processing capacity. In addition, the computational complexity of the proposed algorithm with different numbers of selected batches is given in Fig. 8 (a). We can see that the computation time of the proposed algorithm with 2 11 batches is smaller than the integration time, providing a real-time display. Fig. 8 (b) shows the ISL of the range-Doppler maps obtained by the proposed algorithm with different numbers of batches obtained by 100 random selections and average processing. As can be seen from Fig. 8 (b), the sidelobe level of the proposed algorithm with 2 11 batches is smaller than the batches algorithm with 2 14 batches. Besides, by offsetting the burden of CS by considering 8 times (2 11 vs 2 14 ) fewer cross-correlations for range compression, the proposed algorithm also have comparable processing time with the batches algorithm, as shown in Fig. 8 Finally, about 1 min of the data measured to track the airplane is processed using the proposed algorithm, and the result is shown in Fig. 10. The white rectangle indicates the flying path of the airplane. It can be seen that the landing route of the airplane can be effectively estimated by the proposed algorithm within a 20 dB dynamic range.

B. Large cross-section and slow moving targets: ships

The second demonstration aims at detecting ships entering and leaving the port of Sendai. The objective is to demonstrate the multiple target tracking capability and target separation capacity of the proposed algorithm. Ships are selected because of their lower speed, leading to small Doppler frequencies that are closer to the stationary clutter.

The same analysis sequence as that used above for the airplane measurements is employed for the case of two large ships crossing each other, creating both positive and negative Doppler shifts. Fig. 11 (a) shows the range-Doppler map obtained from the classical batches algorithm. Fig. 11 (b) exhibits the same map obtained using the proposed algorithm, and its result is consistent with that shown in Fig. 11 (a), but the targets are better resolved. Again, background, calculated as the mean value of all batches, is subtracted from each batch, and the Doppler estimation results are obtained from the Fourier transform along the slow time direction. 13 (a) and Fig. 13 (b). To emphasize this result, a short integration time of 64 ms (as opposed to 2.048 s selected previously) was used, decreasing the velocity resolution to a level at which the two targets are mixed when the batches algorithm is used. The 64 ms duration is selected to mix the targets in the Doppler direction, while the practical application would aim at demonstrating separation of targets moving at slightly different speeds, inducing Doppler frequency differences below the inverse of the integration time.

Indeed, despite some suggestions of two maxima in Fig. 14 (blue-dashed curve), the right-most maximum obtained by the classical batches algorithm is not located at the correct speed but at 20 Hz, as opposed to that obtained by the proposed algorithm, which properly identifies the weaker target velocity at 10 Hz, as is shown with the longer integration time in Fig. 12. The extended Doppler frequency range of Fig. 14 is selected to illustrate that the sidelobe peak is stronger than the weaker target when using the classical batches algorithm, as indicated by the black circle in Fig. 14.

The range-Doppler tracking of the two ships during a 2.7min-long acquisition is shown in Fig. 15, where the routes of two ships can be clearly observed. Thus, the new algorithm indeed allows for improved separation of multiple targets characterized by different Doppler shifts.

V. CONCLUSION

We have proposed a novel real-time range-Doppler map generation algorithm for PBR applications. We investigated low-cost PBR consisting of synchronized DVB-T receivers used as general-purpose software-defined radio sources. A low sampling rate of 2.048 Msamples/s is used to demonstrate the performance of the proposed batch CS algorithm for detecting and tracking moving targets. Following the detailed description of the theoretical background and the novel developments, the algorithm was employed using moving airplanes and ships as targets. The proposed algorithm showed an improved range and Doppler sidelobe level reduction compared to the classical batches algorithm while allowing real-time processing, as the computation duration was shorter than the acquisition duration. The improved Doppler resolution helped to improve the target separation capability, as demonstrated in the example of two ships crossing each other. The designed low-cost PBR system and the proposed algorithm have potentials to be applied in airplane navigation and harbor protection. The further study will focus on the extension of the proposed algorithm for the simultaneous estimation of target range, velocity and azimuth information, which will increase the capacity of the PBR system. In addition, CS theory will be employed for range resolution improvement for short-range moving target detection and tracking, which is beneficial for some applications, such as traffic density monitoring.

Correct Doppler identification

  first exponential function in (

2 0 2 =

 22 min || || , . . || || un un surv s t s (15) where ε is the noise level.

  component, mth multipath component, and p-th moving target component in the n-th range profile, respectively. Because the range profiles of DPI and stationary targets show small changes in different batches, the mean value subtraction process is employed to suppress their influence, which is given by 1

where

  NDop is the number of Doppler frequencies. Using the transpose process, we obtain can be considered as a solution obtained by the Fourier transform of the following problem: where (•) T denotes the transpose.

Fig. 1 .

 1 Fig. 1. Flow chart of the proposed range-Doppler map generation algorithm.

  (a) shows the crosscorrelation results of different batches for a total of 2 10 batches. Fig. 3 (b) presents the mean-value subtracted cross-correlation result, with the airplane clearly observed. Fig. 3 indicates that the target can be easily distinguished only after removing DPI and strong stationary target echoes. (a) Cross-correlation results (b) Background subtraction Fig. 3. (a) Cross-correlation results of different batches, and (b) mean values along the slow time is subtracted from all batches, removing static clutter and leaving only the moving target indicated by the black rectangle.

4 .

 4 Range-Doppler map resulting from the Fourier transform along the slow time, providing estimation of the Doppler frequencies introduced by target motion. (a) 2 10 batches and (b) 2 14 batches.

Fig. 5 .

 5 Fig. 5. Sidelobe level of the batches algorithm as a function of the number of batches.

Fig. 6 .

 6 Fig. 6. Proposed algorithm combining the batches algorithm with compressive sensing run on 2 11 randomly selected batches (the white circle indicates target and the white square indicates the outlier).

  (a) ISL against different number of trials (b) Range-Doppler map obtained by averaging over 100 trials Fig. 7. Proposed algorithm combining the batches algorithm with compressive sensing run on 2 11 randomly selected batches: (a) the ISL against different number of trials and (b) range-Doppler map obtained by averaging over 100 random selections (the white circle indicates target).

  (a).

  (a) Computation running time (b) Sidelobe level Fig. 8. (a) Computation running time as a function of selected batch number, demonstrating real-time capability for all cases. (b) Sidelobe level of the proposed algorithm as a function of selected batch number and comparison with the sidelobe level of classical batches algorithm with 2 14 batches. For a more illustrative comparison of the proposed algorithm with the batches algorithm, two cuts are made along the location of the airplane in the Doppler and range directions shown in Fig. 4 (b) and Fig. 7 (b), and the results are shown in Fig. 9. The proposed algorithm has a considerably lower sidelobe level than the batches algorithm, confirming the former's superiority.

  (a) Doppler direction (b) Range direction Fig. 9. Sidelobe level comparison between the batches algorithm and the proposed algorithm in (a) Doppler direction and (b) range direction.

Fig. 10 .

 10 Fig. 10. Airplane tracking in the Doppler-range plane using 1-min data processed by the proposed algorithm.

  (a) Classical batches algorithm (b) The proposed algorithm Fig. 11. Range-Doppler maps obtained by the (a) classical batches algorithm and (b) the proposed algorithm: both ships are visible in a range of 6 km, one exhibiting a radial velocity away from the receiver (stronger target) and the other toward the receiver (weaker target).

Fig. 12

 12 Fig. 12 demonstrates the sidelobe level reduction achieved by the proposed algorithm as a cross-section along the abscissa of Figs. 11 (a) and Fig. 11 (b) for a range in which the two ships cross. Targets moving at different speeds are better separated by the proposed algorithm than by the classical batches algorithm, as shown by the comparison of the two methods in Fig. 13 (a) and Fig. 13 (b).To emphasize this result, a short integration time of 64 ms (as opposed to 2.048 s selected previously) was used, decreasing the velocity resolution to a level at which the two targets are mixed when the batches algorithm is used. The 64 ms duration is selected to mix the targets in the Doppler direction, while the practical application would aim at demonstrating separation of targets moving at slightly different speeds, inducing Doppler frequency differences below the inverse of the integration time.Indeed, despite some suggestions of two maxima in Fig.14(blue-dashed curve), the right-most maximum obtained by the classical batches algorithm is not located at the correct speed but at 20 Hz, as opposed to that obtained by the proposed algorithm, which properly identifies the weaker target velocity at 10 Hz, as is shown with the longer integration time in Fig.12. The extended Doppler frequency range of Fig.14is selected to illustrate that the sidelobe peak is stronger than the weaker target when using the classical batches algorithm, as indicated by the black circle in Fig.14.

Fig. 12 .

 12 Fig. 12. Sidelobe level reduction as observed from cross-sections of the range-Doppler maps for a range in which both ships cross.

Fig. 14 .

 14 Fig. 14. Quantitative analysis of the Doppler resolution improvement and ability to separate targets using the proposed algorithm.

Fig. 15 .

 15 Fig. 15. Tracking of two ships for a duration of 2.7 min: the left ship moves away from the receiver and the right ship moves toward the receiver.

  )

							noise in the surveillance channel.
							The classical way to estimate the bistatic range and Doppler
							frequency of the moving targets is to calculate the CAF between
							the received reference and surveillance signals. The CAF can
							be expressed as
							( , ) f	0 T int	( ) ( ref s t s t surv	)	e	j	2	ft	d t	(3)
	surv s t	surv A s t t P p p a s t t e surv p	j	p f t M m	m a s t t surv n t	m	(2)
	where Asurv is the amplitude of the DPI signal with time delay
	tsurv, which is equal to tref when the reference and surveillance
	antennae are co-located, and M, am, and tm are the number,
	amplitude, and time delay of the multipath echoes, respectively.
	The third term on the right side of (2) denotes the moving target
	components, where fp is the Doppler frequency of the p-th
	moving target; P, ap, and tp are the number, amplitude, and time
	delay of the moving targets; and nsurv(t) denotes the additive

ACKNOWLEDGMENT

This work is supported by JSPS Grant-in-Aid for Scientific Research (A) 26249058. J.-MF is grateful to Tohoku University for granting a position as a three-month visiting scientist to the Center for Northeast Asian Studies. The DVB-T passive radar investigation is partly supported by the French Centre National de la Recherche Scientifique (CNRS) PEPS grant.

This work was supported by JSPS Grant-in-Aid for Scientific Research (A) 26249058. The DVB-T passive radar investigation is partly supported by the French Centre National de la Recherche Scientifique (CNRS) PEPS grant.


Stavros Meskos
email: stameskos@imbm.bas.bg

Stefan Stefanov
email: stefanov@imbm.bas.bgc

Dimitris Valougeorgis

DSMC Calculations of Binary Gas Mixing in Simple Micro-Sized Configurations

The Direct Simulation Monte Carlo (DSMC) method is employed to model and simulate the mixing of two non-reactive gas streams in a simple twodimensional micro-mixer. The efficiency of the mixing process in a mixer is related to the mixing length, which is the distance from the inlet where the mixture is found to be homogeneous. A new method to measure homogeneity is presented. It has been found that the boundary conditions at the walls and the inlets, the molecular models and the physical properties of the gas species may influence the mixing process at different scales. Their level of influence on the efficiency of the mixer is presented.

INTRODUCTION

Gas mixing in micro-scale is a problem of theoretical and industrial interest for development and optimal design of gaseous MEMS/NEMS (Micro/Nano-Electro-Mechanical Systems) devices [1]. The study of the effects on the mixing process for various configurations and different boundary conditions is essential for the design of these devices. The computational analysis of gaseous flows in MEMS devices operating under different thermal conditions in non-equilibrium flow regimes cannot be based on classical continuum models of fluid motion because the continuum assumption that the flow is locally in the near-equilibrium state is no longer valid. Thus, numerical methods that solve the Boltzmann equation such as the Discrete Velocity Method and the Lattice Boltzmann Method, as well as particle based methods such as the Molecular Dynamics and the Direct Simulation Monte Carlo methods are used instead. The Direct Simulation Monte Carlo (DSMC) method [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] is arguably the most common and is employed in this work.

One of the first investigations of gas mixing in micro-configurations has been performed by Yan & Farouk [START_REF] Yan | [END_REF]. They investigated the dependence of the mixing length on the pressure ratio between the inlets and the outlets, between the inlets alone and on the inlet velocities. The tested configuration consisted of two parallel gas streams (H 2 -O 2 ) entering the mixing chamber. In the inlet region the two streams are separated by a splitter plate. It has been found that the mixing length is increased by increasing the pressure ratio between inlets and outlets and between the two inlets. In all simulations the outlet pressure is kept at 50 kPa. The same flow setup has been also examined by Wang & Li [4] testing the mixing of CO and N 2 while the outlet pressure is kept at zero (no backflow). They introduced the terms of relative density difference and mixing coefficient in order to describe better the mixing process. It has been found that the mixing length is inversely proportional to the gas temperature and the Knudsen number, while it is proportional to the Mach number. They have also shown that the wall characteristics have little effect on the mixing length when the rest of the properties are the same. Hassan & Le [5] studied the mixing of CO and N 2 at a T-shape micro-mixer by also keeping the outlet pressure low enough to prevent back-flow. It has been found that at higher Knudsen number the mixing length is reduced, while increasing the inlet pressure resulted in an increased mixing length. Also, increasing the flow rate and the wall temperatures reduced the mixing length, with the wall temperature effect to be more significant. Most recently, Darbandi and Sabouri [6] have investigated CO and N 2 mixing at an angle by applying Y-shaped inlets. They have also studied the effect of replacing a "larger" micro-mixer with many "smaller" ones.

In this work, following [4], the mixing of two parallel gas streams entering the microchannel is simulated using the DSMC method. Initially, the influence of different parameters, such as temperature and accommodation coefficients at walls, inlet pressures and molecular models, on the mixing process is investigated using CO and N 2 which have the same molecular mass. Then, in order to investigate the mixing process of gases with different molecular masses the mixing of the gas pairs of He-Xe, Ne-Ar, He-Ar and He-Ne is considered.

PROBLEM DESCRIPTION

The basic configuration of the micro-mixer with dimensions L/H = 10, H/d = 1, H = 1μm is shown in Fig. 1. The two gases enter from the left upper and lower inlets respectively. As "one" is considered the lighter species, which enters from the upper inlet, unless stated otherwise. For the case of CO and N 2 mixing the former gas is considered as one. Pressure at each inlet is 0.2 atm and temperature is set at 300 K, while vacuum is assumed at the outlet. Initially, vacuum conditions are assumed also inside the mixer. The accommodation coefficients at all walls are equal to unity, (fully diffuse reflection) and the inlet bulk velocity is zero. The area after the inlets is the "mixing chamber". The two species start mixing when they enter the mixing chamber. The mixing length is defined by the distance from the beginning of mixing chamber to a distance where the two species are found to be fully mixed. The mixture is considered as fully mixed when the concentration of the species across the channel's height is equal to the ratio between the upper and lower fluxes. One way to calculate the mixing length is by measuring the homogeneity of the mixture at different lengths across the channel [4].The proposed method reads as follows:

,, ,, max( , )

j up j low aa j a j up j low aa nn n    (1)
Here, j a  is the relative density difference, which will be referred to as "RDD" in this paper, of species a at column j . Also, , j up a n and , j low a n correspond to the number densities of the same species at the same column as measured at the two cells adjacent to the upper and lower walls of the main channel (L) respectively. Across the length of the channel it is expected that the RDD of each species will reach a minimum value. In an ideal case it would be equal to zero but for practical reasons in this work it is set at 0.5%. The length at which each species is found to reach this minimum value it is the point that each species is considered to be uniformly distributed across that column. The lengths that this happens for each species might coincide or not. The larger of those two is chosen as mixing length.

Based on the same concept an alternative method is proposed in this work, where all intermediate cells are considered:

, , , 1

1 j j avg N a i a j a j avg i a nn Nn      (2) , , 1 1 N j avg j a a i i nn N    (3)
Here, 1,.., iN 

is the number of cells at column j and , j avg a n is the average number density over the cells of column j . This approach is considered as more accurate since the average of the deviations of each cell number density of the column average is calculated.

In Fig. 2, the evolution of the RDD over the length of the mixing chamber computed with the present method and the one proposed in [4] is plotted. It is observed that the RDD evolution of both species based on the two methods is similar and this is justified by the same molecular masses of CO and N 2 . However, with the Wang & Li method the statistical noise of DSMC is very significant, since the computation is based only on the extreme cells across the channel. For both cases the minimum value of RDD is about 0.005 or 0.5% at wall length equal to 4.8 μm or at a mixing length equal to 2.8 μm. The sudden increase of the slope at the end of the RDD evolution (obtained with the proposed method) indicates enhanced separation, which is justified by the imposed vacuum at the outlet.

COMPUTATIONAL CONSIDARATIONS

The DSMC method, which is used in the present study, is a well-known particle based stochastic method that simulates gas flows very efficiently in transition and rarefied regimes. The main principle of the method is the decupling of the collisions and movement of particles. This is achieved by considering a time step smaller than the time a particle travels a mean free path at most probable velocity. Each simulated particle represents a huge amount of real molecules, usually greater than 10 18 . The computational domain is divided into cells. Particles are moving through the cells, collisions are occurring only between particles within the same cells and finally molecular speeds and number densities are sampled at each cell. This process is repeated at each time step. The simulation continues in time until a steady state is reached. There are two possibilities of averaging. The first is done by time averaging a large number of time steps after reaching the steady state. This averaging method is sufficient when only the steady state solution is of interest. However, in order to derive accurate results for the transient period, an ensemble averaging method is required as well, where many independent simulations are performed and results are additionally averaged at the same time steps by ensemble averaging over all simulations. An in-house DSMC code has been developed and employed in this work. For a rigorous mathematical description of the method, the reader is referred to [7].In Table 1, the basic DSMC configuration that has been applied in all simulations, unless stated otherwise, is presented.

RESULTS AND DISCUSSION

The influence of the wall accommodation coefficient on the mixing length is depicted at Figures 3 (a) and (b). Figure 3(a) corresponds to the case where the accommodation coefficient is changed at all walls, while in Fig. 3(b)the accommodation coefficient is altered only at the upper wall while in all other walls is fixed to one. Results are provided with minimum relative density difference of 0.5% and 1% marked in red and blue, respectively. In Fig. 3(a) it is seen that the mixing length is greatly reduced when moving from fully specular to fully diffuse walls. In addition, the requirement for a higher accuracy of a 0.5% results in a larger mixing length (almost up to 1 μm). This becomes more evident at Fig 3(b), where when the RDD is 1%, the mixing does not reach the required threshold value for accommodation coefficient values in the range below 0.3, while when the RDD is 0.5% this is happening in a much wider range of values less than 0.7. It is evident that the wall material and the type of gas-surface interactions have a significant influence on the mixing process. Table 2 compares the mixing length for different temperatures at walls, while the temperature at the inlets is fixed at 300 K. While the temperature is increased at all walls at 1.5 and 2 times, from 300K to 450K and 600K respectively, only a very small decrease in the mixing length is observed. On the third column the optimum RDD achieved at this length is shown and its analysis suggests that although the length is reduced the relative density difference is increased which makes the mixture less homogenous. This remark is further strengthened by the average RDD at the end of the mixing chamber tabulated in the last column of Table 1. As the wall temperature is increased the separation effect of the mixture becomes stronger. [8]. In addition, the Variable Hard Sphere (VHS) and the Hard Sphere (HS) models, proposed by Bird [START_REF] Bird | MONTE-CARLO SIMULATION IN AN ENGINEERING CONTEXT[END_REF] have also been investigated and the results have been compared with the ones obtained by the VSS model. The readers are referred to [START_REF] Bird | The DSMC Method Version 1.2[END_REF] where a detailed and intuitive description of the phenomenological models is provided. Table 3 depicts the mixing length variations over the different models. It is interesting that even for the flow setups with isothermal walls there is a significant difference between the models. In order to better understand the origin of this difference a comparison of the corresponding velocity profiles is conducted. The velocity profiles of the mixture for the VSS case are plotted at the beginning, the middle and the end of the mixing chamber in Fig. 4 as the reference case. Then, a comparison is done at three different points on the center line at each length for each model respectively and Table 4 tabulates the differences in the velocities.

FIGURE 4.

Velocity profiles of the mixture, for the VSS case, at the start, middle and end of the mixing chamber respectively. Velocity reference is the most probable velocity of CO at 300 K (422.08 m/s). Although the differences are relatively small, the influence on the mixing process is large, which means that the velocity of the flow field is also a significant parameter. The authors in [START_REF] Yan | [END_REF]4] arrived in the same conclusion by investigating different inlet bulk velocities. This result supports the Koura and Matsumoto [8] conclusion that the VSS should be preferred in cases of gas mixture flows as a more reliable model.

The effect of the molecular mass difference has also been investigated by considering species with different molecular masses. In Table 5 the results of mixing lengths and RDDs are provided for the following combinations: He-Xe, Ne-Ar, He-Ar and He-Ne. Compared to the base case where the two species with similar masses are mixed, now the relative density difference evolves differently. Several interesting phenomena take place. The lighter species can achieve a much better uniform distribution having RDDs 0.5%-0.7% at a length of about 4.3 μm, while the heavier species never reach that good distribution, with a minimum RDD of about 1% and moving up to 2.4%. In addition, the lighter species is always influenced by the heavier one loosing very quickly its uniformity. Furthermore, the final mixing length is similar in all the cases. The main difference and influence of the heavier species lies on the RDD. This means that different cases may result in the same mixing lengths but with different homogeneity at those corresponding lengths. It is important to note that these species will never be "fully mixed" and only a best possible mixing can be achieved, depending on the boundary conditions. This can be better understood from Fig. 5, where the evolution of RDD for the case of He-Xe, are given with colors (red for He and blue for Xe). The two species are distributed at different rate and their RDD evolution lines never cross or coincide. In addition, at the same figure the RDD evolution of He-Xe mixture is plotted with green for He and black for Xe in a mixing chamber where its length has been increased by 4 μm (total length of the mixer becomes 12 μm). This is performed in order to investigate the influence of the length on the final results. The behavior is identical. In both cases both species reach the minimum value of RDD at the same point and after that they expand towards vacuum in a similar manner.

CONCLUSIONS

The mixing process of two parallel gas streams of monoatomic, non-reactive species have been investigated in terms of the mixing length and homogeneity of the mixture. The boundary conditions at the walls and the inlets of the mixer, the molecular model and the molecular masses of the species have been shown to have different degree of influence to the mixing process. In order to describe correctly that process and compare the efficiency of different micro-mixers it is necessary to define along with the mixing length the corresponding measure of homogeneity. In this work a new type of measurement is proposed by evaluating how uniformly each species is distributed in space. By accepting a very small value, usually about 0.5% the condition of "fully mixed" flow practically is considered as fulfilled. Even a small change in that value can lead to errors. More predominantly for the realistic case of mixing species with different masses this value does not stay fixed; rather it varies differently for each species, giving as well a good indication of the separation effects that take place. Moving from specular to diffuse walls the mixing process is greatly improved. Increasing the wall temperature, with fixed inlet gas temperature, results in a small decrease of the mixing length on the account of the homogeneity. Simulations performed using three different intermolecular models of VSS, VHS and HS resulted in differences in the mixing length that should not be neglected. An explanation for that is related to the observation that the application of different molecular interaction models reflects in small differences in velocity fields of the gas components. Finally, the difference between molecular masses of the two species influences the mixing process mainly on the mixture homogeneity, while the mixing length remains almost unaffected.

FIGURE 1 .

 1 FIGURE 1. Schematic of the parallel gas micro-mixer.

FIGURE 2 .

 2 FIGURE 2. Evolution of relative density difference of CO and N 2 calculated with both methods. At the upper right the evolutions based on the two cases are zoomed.

FIGURE 3 .

 3 FIGURE 3. Mixing length variation over different accommodation coefficients  at walls: (a) all walls have the same

FIGURE 5 .

 5 FIGURE 5. Evolution of relative density difference for the case of He-Xe for mixing chambers with total length of 8 μm (red and blue) and 12 μm (green and black).

TABLE 1 .

 1 DSMC properties

	Property	Value
	Collision scheme	NTC
	Molecular model	VSS (monoatomic)
	Time averaging	20 kinetic steps
	Total samples	500 (1 sample = 20 kin. Steps)
	Ensemble averaging	50 (simulations)

TABLE 2 .

 2 Influence of the wall temperature on the mixing length.

	Temperature [K]	Mixing length [μm]	RDD	Average RDD
	300	2.8	0.57%	0.76%
	450	2.7	0.67%	0.98%
	600	2.5	0.69%	1%
	All simulations have been performed by employing the Variable Soft Sphere (VSS) molecular model, proposed
	by Koura and Matsumoto

TABLE 3 .

 3 Influence of the molecular model on the mixing length.

	Molecular	Mixing length	Relative Error
	Model	[μm]	VSS/VHS VSS/HS
	HS	3.3
	VHS	3.15	12.4%	17.4%
	VSS	2.8

TABLE 4 .

 4 Comparison of velocities for different physical models.

	Length in	Velocity (  422.08 m/s)		Relative Error
	Mixing
	chamber
	[μm]	VSS	VHS	HS	VSS/VHS VSS/HS VHS/HS
	2	0.1837 0.1885 0.1858	2.61%	1.14%	1.43%
	5	0.3133 0.3232 0.3172	3.16%	1.24%	1.88%
	8	0.7692 0.7785 0.7682	1.21%	0.13%	1.32%

TABLE 5 .

 5 Mixing lengths and RDDs for species with different physical properties

	Species	RDD1	Point 1 [μm]	RDD2	Point 1 [μm]	Mixing length [μm]	m1 / m2 d1 / d2
	He-Xe	0.6% -1%	4.23	1.4% -2.4%	4.16	2.23	0.03	0.404
	Ne-Ar	0.7% -1.1%	4.23	1%	4.16	2.23	0.506	0.659
	He-Ar	0.5% -1.1%	4.3	1.1% -1.6%	2.56	2.56	0.1	0.559
	He-Ne	0.7% -1%	4.3	1% -1.5%	2.5	2.5	0.2	0.845

ACKNOWLEDGMENTS

This project has received funding from the European Union's Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Sklodowska-Curie Grant Agreement No. 643095.


Yun Chen

Azure Hansen

Rodolphe Boudot

Moshe Shuker

Eugene Ivanov

Gregory Hoth

Bruno Pelle

John Kitching

Elizabeth A Donley

Two-Dimensional

Two-dimensional rotation measurement with point source interferometry

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Two-Dimensional Rotation Measurement with Point Source Atom Interferometry

Yun-Jhih Chen 1,2 , Azure Hansen 1 , Rodolphe Boudot 1,3 , Moshe Shuker 1,2 , Eugene Ivanov 1,4 , Gregory W. Hoth Abstract-In analogy with a spinning rotor gyroscope, which precesses in response to a torque in the plane transverse to the axis of spin, point source atom interferometry (PSI) senses the components of rotation in the plane transverse to the direction of a pair of counter-propagating Raman laser beams. Based on the PSI technique, we demonstrate the two-dimensional measurement of a rotation vector projected into a plane. We characterize the sensitivity of the measurement, including the magnitude and the direction of the rotation vector component.

Keywords-gyroscope; atom interferometer; quantum sensor

Rotation measurements with light pulse atom interferometry are typically more complex than acceleration measurements because both acceleration and rotation contribute to the interferometer phase shift. Two counter-propagating cold-atom sources or more complicated atom-optics are often used to obtain the rotation measurement from the interferometer phase shift [START_REF] Canuel | Six-axis inertial sensor using cold-atom interferometry[END_REF]. In contrast, point source atom interferometry (PSI) [START_REF] Dickerson | Multiaxis inertial sensing with long-time point source atom interferometry[END_REF] makes use of the thermal velocity spread of an isotropically expanding cloud of cold atoms to construct the rotation measurement. While dramatically simplifying the architecture of an atom interferometer gyroscope, the PSI's unique capability to simultaneously measure two components of the rotation vector with moderate to high bandwidth also opens up new applications. For instance, it may permit measuring timevarying rotation vectors, such as precessions, with atom interferometers.

In PSI, the most common Raman beamsplitter-mirrorbeamsplitter pulse sequence is applied to a cold-atom cloud which expands within a pair of counter-propagating Raman laser beams. If the final cloud is much larger than the initial cloud, the latter is an approximation to a "point source" and the position of an atom in the final cloud is related to the velocity of the atom by 𝑟 ⃑ = 𝑣 ⃑𝑇 ex , where 𝑇 ex is the total expansion time of the cloud. This position-velocity correlation maps the interferometer phase shift of each atom into space as Φ(𝑟 ⃑) = 𝑘 ⃑⃑ eff • 𝑎 ⃑𝑇 2 + 𝑘 ⃑⃑ Ω • 𝑟 ⃑, where the first term is the acceleration phase and the second term is due to a phase gradient induced by rotation. The phase gradient is 𝑘 ⃑⃑ Ω = (𝑘 ⃑⃑ eff × Ω ⃑⃑⃑ )2𝑇 2 /𝑇 ex , where Ω ⃑⃑⃑ is the angular velocity, 𝑘 ⃑⃑ eff is the effective Raman wavevector, and 𝑇 is the time between the Raman pulses. The population ratio is a sinusoidal function of the interferometer phase, so the phase gradient creates a spatial fringe pattern in the final cloud (Figure 1). The projection of Ω ⃑⃑⃑ into the plane transverse to 𝑘 ⃑⃑ eff is constructed from the fringes. The number of the fringes is proportional to the magnitude of the rotation projection. The orientation of the fringes is in the direction of the rotation projection. The acceleration in the direction of 𝑘 ⃑⃑ eff translates the fringes without affecting the number and orientation of the fringes. With a single cloud of laser cooled 87 Rb atoms, we have previously demonstrated rotation measurements with PSI in the regime where the final cloud is only a few times larger than the initial cloud. We have characterized a systematic error due to this non-ideal point source [START_REF] Hoth | Point source atom interferometry with a cloud of nite size[END_REF]. In the current work, we focus on the two-dimensional rotation measurement of PSI. We demonstrate the measurement of a rotation vector projected onto the plane normal to the Raman laser beams and characterize the sensitivities of the measurement, including the magnitude and direction of the rotation vector projection in the plane (Figure 2).

In the presentation, we will present the details of the 2D rotation measurement with PSI and discuss our ongoing work on investigating the technical performance limitations and systematic errors in our compact setup, in which we have such a non-ideal point source. We will also discuss the applications of the PSI technique; for example, the 2D rotation measurement with PSI can be used to find geographic north (Figure 3). ⁄ power fits from 𝜏 = 1 to 10 s. For a 1 s averaging time, the sensitivity of the magnitude is 0.033 °/s and that of the direction is 0.27 °. We estimate that, with an ideal point source limited only by atom shot noise, for 1 s the ultimate sensitivities could be as low as 1 × 10 -4 °/s and 3 × 10 -3 °, respectively.

Fig. 1 .

 1 Fig. 1. Rotation fringes obtained with a simulated rotation. Images are autoscaled with red indicating population maxima and blue population minima. The Raman laser beams propagate parallel to the z-axis and the cold-atom cloud is imaged in the xy-plane. Here, 𝑇 = 8 ms and 𝑇 ex = 26 ms.

Fig. 2 .

 2 Fig. 2. Allan deviation of the magnitude and direction of the rotation vector measurement in a plane with simulated rotation. Here, 𝑇 = 8 ms and 𝑇 ex = 26 ms. Black lines are 𝜏 -1 2⁄ power fits from 𝜏 = 1 to 10 s. For a 1 s averaging time, the sensitivity of the magnitude is 0.033 °/s and that of the direction is 0.27 °. We estimate that, with an ideal point source limited only by atom shot noise, for 1 s the ultimate sensitivities could be as low as 1 × 10 -4 °/s and 3 × 10 -3 °, respectively.

Fig. 3 .

 3 Fig. 3. Illustration of gyrocompassing with PSI. When the direction of the Raman laser beams (𝑘 ⃑⃑ eff ) is normal to the surface of Earth, PSI senses the projection (Ω ⊥ ) of the Earth rotation (Ω ⃑⃑⃑ ) onto the Earth's surface. Since Ω ⊥ points to the geographic north, the rotation fringes due to Earth's rotation points to the geographic north.

  1,2* , Bruno Pelle 1,2 † , John Kitching 1 , and Elizabeth A. Donley 1 1. Time and Frequency Division, National Institute of Standards and Technology, Boulder, CO 80305 USA 2. University of Colorado Boulder, Boulder, CO 80309 USA 3. FEMTO-ST, CNRS, 26 Rue de l'Epitaphe, 25030, Besançon, France 4. University of Western Australia, Perth, WA, Australia

ACKNOWLEDGMENT

We thank Mark Kasevich for helpful discussions. A. H. was supported under an NRC Research Associateship award at NIST. R. B. was supported by the NIST Guest Researcher Program and the Délégation Générale de l'Armement (DGA). This work was funded by NIST, a U. S. government agency, and is not subject to copyright.


Moshe Shuker

James Wesley Pollock

John Kitching

Elizabeth A Donley

Rodolphe Boudot

Valeriy I Yudin

Alexey V Taichenachev

Light-Shift Suppression with Novel Variants of Adaptive Ramsey Spectroscopy

Keywords: atomic clocks, light shifts, Ramsey spectroscopy I

come   L'archive ouverte pluridisciplinaire

BACKGROUND

Ramsey's method of separated oscillating fields [START_REF] Ramsey | A Molecular Beam Resonance Method with Separated Oscillating Fields[END_REF] is a powerful technique that is widely used for the measurement of atomic transition frequencies. With standard Ramsey spectroscopy, two short interrogation pulses are separated by a dark period, with the first pulse generating the atomic coherence and the second pulse probing its phase. Although the atoms spend most of their time in the dark, non-negligible residual sensitivity to light shifts is observed due to shifts that arise from the interrogation pulses. These light shifts cause systematic frequency shifts and instability in optical atomic clocks [START_REF] Ludlow | Optical atomic clocks[END_REF] as well as clocks based on coherent population trapping (CPT), which probe microwave transitions optically [START_REF] Arimondo | Coherent population trapping in laser spectroscopy[END_REF][START_REF] Vanier | Atomic clocks based on coherent population trapping: a review[END_REF][START_REF] Shah | Advances in Coherent Population Trapping Atomic Clocks[END_REF].

Sophisticated interrogation protocols based on Ramsey's method and tailored composite laser pulses have been proposed and demonstrated that can eliminate probe-induced frequency shifts (for a review, see [START_REF] Zanon-Willette | Composite laser-pulses spectroscopy for high-accuracy optical clocks: a review of recent progress and perspectives[END_REF]). Key to these techniques is that the light shift magnitude for standard Ramsey spectroscopy depends on the duration of the Ramsey dark period, scaling as 1 / 𝑇𝑇. The light shift also depends on several parameters related to the pulses themselves, including the interrogation frequency, phase, and/or pulse length. With a technique called hyper-Ramsey spectroscopy [START_REF] Yudin | Hyper-Ramsey spectroscopy of optical clock transitions[END_REF], the pulse parameters can be tuned to form composite pulses that compensate for the light shifts. One can determine the appropriate form of the composite pulses by performing Ramsey frequency measurements with various composite pulses for different values of 𝑇𝑇, and tuning the pulse parameters until the central fringe frequency becomes independent of 𝑇𝑇, at which point light shift vanishes.

Autobalanced Ramsey (ABR) spectroscopy is a technique that takes the removal of light shifts one step further by actively adjusting the composite pulses to stabilize against the shifts [START_REF] Sanner | Autobalanced Ramsey Spectroscopy[END_REF]. With ABR, two error signals derived from interleaved Ramsey sequences of different dark periods drive two feedback loops, with one feedback loop applied to the clock frequency and one feedback loop used to control a phase correction applied between the pulses. Autobalanced Ramsey spectroscopy was proposed and originally demonstrated in a ytterbium single-ion clock [START_REF] Sanner | Autobalanced Ramsey Spectroscopy[END_REF], in which a phase correction was applied to the second Ramsey pulse, thus cancelling the light shift. The method was later applied to vapor-cell clocks based on CPT to show a significant improvement in mid-term stability [START_REF] Hafiz | Toward a High-Stability Coherent Population Trapping Cs Vapor-Cell Atomic Clock Using Autobalanced Ramsey Spectroscopy[END_REF][START_REF] Hafiz | Symmetric autobalanced Ramsey interrogation for high-performance vapor-cell atomic clock[END_REF].

The ABR technique has also been generalized theoretically to show that a variety of secondary variables (or concomitant parameters) in addition to a controlled phase correction can be used for the composite pulses, including external frequencystep compensation and variable second-pulse duration [START_REF] Yudin | Generalized Autobalanced Ramsey Spectroscopy of Clock Transitions[END_REF].

II. DISPLACED FREQUENCY JUMP RAMSEY

A new technique called displaced frequency-jump Ramsey spectroscopy (DFJR) was recently demonstrated [START_REF] Shuker | Ramsey Spectroscopy with Displaced Frequency Jumps[END_REF], which is related to but different from ABR in a few important ways. The fringes in Figure 1 illustrate the concept of the DFJR technique.

Like ABR, it uses two interleaved Ramsey sequences with different dark periods and uses two control loops, but it does not use adjusted composite pulses in one of the control loops. Instead, the controlled concomitant parameter, 𝛼𝛼, is a scaling factor for a finely adjusted frequency displacement that is inversely proportional to the dark period and applied to the frequency jumps that probe the sides of the central Ramsey fringe: 𝛥𝛥𝑓𝑓 𝑇𝑇 = 𝛼𝛼 ⁄ 𝑇𝑇 . Thus, the probe frequency of the local oscillator (LO) modulates about the central fringe in the Ramsey spectrum like usual in a Ramsey-style clock, but with a servoed frequency displacement having a 1 / 𝑇𝑇 dependence determined from interleaved sequences of different dark periods. This displacement is tracked in software, and when it is subtracted from central Ramsey fringe frequency, the result is the unshifted clock frequency, 𝑓𝑓 𝑐𝑐 . The value of the light shift is measured by the scaled frequency displacement 𝛥𝛥𝑓𝑓 𝑇𝑇 . The DFJR method has been shown to suppress light shifts in the cold-atom CPT clock developed at NIST [START_REF] Esnault | Cold-atom double-Λ coherent population trapping clock[END_REF][START_REF] Liu | Low-Drift Coherent Population Trapping Clock Based on Laser-Cooled Atoms and High-Coherence Excitation Fields[END_REF][START_REF] Pollock | ac Stark shifts of dark resonances probed with Ramsey spectroscopy[END_REF] by more than one order of magnitude [START_REF] Shuker | Ramsey Spectroscopy with Displaced Frequency Jumps[END_REF] to the 1 × 10 -12 level or lower in our experiment. There are a few technical factors that make DFJR attractive in practice. The control parameters, 𝛼𝛼 and 𝑓𝑓 𝑐𝑐 , are controlled with a single physical variable -the LO frequency. Modulation of the LO phase is not needed, which allows for simple implementation without systematics that can arise from modulating the oscillator phase. Also, unlike one variant of generalized ABR [START_REF] Yudin | Generalized Autobalanced Ramsey Spectroscopy of Clock Transitions[END_REF] where frequency steps are applied during the Ramsey pulses, here the LO frequency remains constant during the Ramsey cycle, and no rapid control of the LO frequency is required.

III. COMBINED ERROR SIGNAL RAMSEY

It has also been shown theoretically [START_REF] Yudin | Combined error signal in Ramsey spectroscopy of clock transitions[END_REF] and experimentally [START_REF] Shuker | Reduction of light shifts in Ramsey spectroscopy with a combined error signal[END_REF] that probe-induced light shifts can be removed by combining two error signals into a single error signal with a method called combined error signal Ramsey spectroscopy (CES) [START_REF] Yudin | Combined error signal in Ramsey spectroscopy of clock transitions[END_REF]. The error signal steers the clock frequency to its unshifted value without requiring a concomitant parameter. Like before, the two error signals are generated with alternating Ramsey sequences. To combine the signals with the appropriate weights, the error signals are normalized by a factor that equalizes the Ramsey fringe amplitude and then subtracted. The method is simple since there is a single control loop with the clock frequency as a single control parameter.

A challenge with implementing this method is to determine the appropriate normalization factors, but there is a natural way to do this with NIST's cold-atom CPT clock, in which the main decay mechanism arises from atoms escaping the probe volume after they are released from the optical molasses. So, the normalization process in this case can easily be accomplished by inserting a delay before the shorter Ramsey pulse in the sequence such that the probe pulses (the 2 nd Ramsey pulses) occur after the same release period from the molasses for the long and short Ramsey sequences. The CES method has also been shown to suppress light shifts by at least one order of magnitude over traditional Ramsey-CPT spectroscopy [START_REF] Shuker | Reduction of light shifts in Ramsey spectroscopy with a combined error signal[END_REF].

IV. OUTLOOK

So far, our work has focused on quantifying light-shift suppression in our cold-atom clock at NIST. An ongoing important focus of research will be to measure how much the long-term stability of CPT clocks can be improved by these techniques, and the apparatus is currently undergoing an upgrade to be able to perform these measurements.

Autobalanced Ramsey and other adaptive spectroscopy methods are a new field with many interesting potential applications. In principle, it should be possible to remove systematics other than light shifts when there is a model that closely predicts the scaling of the shift versus 𝑇𝑇.

Fig. 1 .

 1 Fig. 1. Simulated light-shifted central Ramsey fringes for two different dark periods illustrating the concept of the DFJR technique. In the top figure, the Ramsey period was 16 ms, giving a fringe width at half maximum of 31.25 Hz. The light shift, 𝛥𝛥𝑓𝑓 𝑇𝑇1 = 𝛼𝛼 ⁄ 𝑇𝑇 1 , was assumed to be 7 Hz (made large for illustrative purposes). In the bottom figure, the Ramsey period was 4 ms, the fringe width was 125 Hz, and since the Ramsey period was four times shorter than in the top figure, the light shift, 28 Hz, was four times larger. The frequency jumps for the servo modulation, fJ,1(2) = 1/(4⸱T1(2)), are different for the two dark periods and are offset with the second control loop by the frequency displacement ΔfT. ΔfT is found along with the clock frequency fc by interleaving Ramsey sequences with two different dark periods. In effect, the short Ramsey sequences steer the frequency displacement 𝛥𝛥𝑓𝑓 𝑇𝑇 = 𝛼𝛼 ⁄ 𝑇𝑇 to compensate for the light shift and the long Ramsey sequences steer the clock frequency fc to the unshifted frequency f0.

ACKNOWLEDGMENT

The authors thank V. Maurice, L. Stern, and M. Lombardi for helpful suggestions on this manuscript. V. Maurice is also gratefully acknowledged for writing the software used to perform the experiments. V. I. Yudin was supported by the Russian Science Foundation (No. 16-12-10147), Ministry of Education and Science of the Russian Federation (No. 3.1326.2017/4.6), and Russian Foundation for Basic Research (No. 17-02-00570). A. V. Taichenachev was supported from Russian Science Foundation (No. 16-12-00054) and Russian Foundation for Basic Research (No. 18-02-00822). R. Boudot was supported by the NIST Guest Researcher Fellowship and Délégation Générale de l'Armement (DGA).