Wassim Khaddour

Foudil Dadouche

Wilfried Uhring
email: wilfried.uhring@unistra.fr

Vincent Frick

Morgan Madec

Design Methodology and Timing Considerations for implementing a TDC on a Cyclone V FPGA Target

Keywords: Time-to-Digital Converter (TDC), Field-Programmable Gate Array (FPGA), Carry Chain Logic, routing delay

There are hundreds of research publications that theoretically discuss the implementation of Tapped Delay Line based Time to Digital Converters (TDL TDCs) on Field-Programmable Gate Array (FPGA) targets. However, most of these works do not cover the timing issues that will be encountered mostly due to the routing delays. The purpose of this work is to highlight the main timing issues that should be considered when implementing TDCs in FPGA targets and propose practical approaches to overcome these issues. As a study case, a full design methodology of a TDC on a Cyclone V FPGA target is presented in this work.

I. INTRODUCTION

Time to Digital Converters (TDC) are widely used in numerous applications that require the measurement of precise time intervals, such as the time of flight (TOF) based applications, the time-resolved fluorescence measurement, LIDAR applications, fluorescence lifetime imaging, 3D active imaging, time correlated photon counting, etc. The role of the TDC in these applications is to measure the time of flight of the photons, i.e. the time interval that separates the emission of a light pulse and the detection of a photon. There are many architectures of TDCs such as: The Tapped-Delay Line (TDL), Vernier Delay Line (VDL), Delay Locked Loop (DLL) and the phased clocks. These architectures are usually designed as Application-Specific Integrated Circuits (ASIC), which are high-performance solutions, but at high measurement rate, the data transfer become difficult to handle. Moreover, the downstream processing could be more efficiently integrated in Field-Programmable Gate Arrays (FPGAs), thus implementing the TDC's within the FPGA is a better solution [START_REF] Sui | A 2.3-ps RMS resolution time-to-digital converter implemented in a low-cost Cyclone V FPGA[END_REF].

Over the last years, many research works covered the different architectures for implementing TDCs in FPGA. Most of these works focused on the improvements that could be made to the TDCs in terms of resolution and linearity [START_REF] Sui | A 2.3-ps RMS resolution time-to-digital converter implemented in a low-cost Cyclone V FPGA[END_REF]- [START_REF] Machado | Recent developments and challenges in FPGA-based time-to-digital converters[END_REF], but they did not cover the different timing issues that should be taken in consideration when implementing a TDC on FPGA. This work aims to introduce a full design methodology of a TDC on a low cost Intel FPGA Cyclone V development kit, covering the most important timing considerations when integrating this TDC in a global system.

The general functional principle of the TDC is covered in section II, as well as the TDL architecture used to build the studied TDC. Section III addresses the principal key points for implementing a TDC in a Cyclone V target. In section IV, a full design methodology of a TDL is provided. This section presents the global architecture of the studied design and highlights the most critical timing issues that should be considered. Finally, section V presents the conclusion.

II. TDC'S OPERATING PRINCIPLE

A. TDC's general operating principle

TDCs are widely developed to measure precise time intervals between two events represented by two signals, named START and STOP signals, and to convert these intervals into binary values. Most of the TDCs structures consist of three parts; a coarse block and two fine blocks. The coarse block is a simple counter working at the system clock. It counts the number of the entire clock periods between the rising edge of the START and the STOP signals. The role of the fine blocks is to measure the two accurate time intervals; (Tfine1) between the rising edge of the START signal and the following clock's rising edge, and (Tfine2) between the rising edge of the STOP signal and the following clock's rising edge. As illustrated in the timing diagram presented in Fig. 1, the measured time (Tmeasured) is composed of three intervals; Tcoarse, which is equal to the coarse counter multiplied by the system clock's period, Tfine1 and Tfine2. This timing diagram shows that the measured time is given by the following equation [START_REF] Dadouche | Design methodology of TDC on low cost FPGA targets[END_REF].

Tmeasured = Tfine1 + Tcoarse -Tfine2 (1)
= Tfine1 + N.Tclk -Tfine2

This coarse-fine architecture increases the dynamic range of the TDC and improves its resolution.

B. TDC's fine block structure

The fine block is the core of the TDC. The TDC studied in this work is based on the TDL structure. In this architecture, the hit signal, which is the START or the STOP signal, propagates through delay elements, the output of these elements is sampled at the system global clock and encoded to a fine timestamp. The number of delay elements is defined by the clock period and the propagation time of an element. It is common to use the carry propagation-chains included in FPGA platforms to build a TDL TDC structure. The interest of using these chains is that their propagation time is much shorter than the propagation time of the general routing networks which improves the TDC's resolution [START_REF] Dadouche | Design methodology of TDC on low cost FPGA targets[END_REF] [START_REF] Wu | Firmware-only implementation of timetodigital converter (TDC) in field-programmable gate array (FPGA)[END_REF].

As explained in [START_REF] Dadouche | Design methodology of TDC on low cost FPGA targets[END_REF] and [START_REF] Dadouche | Design implementation and characterization of time-to-digital converter on low-cost FPGA[END_REF], to take advantage of the carry propagation-chains of low cost Intel FPGA targets, it is possible to make use of adders. In such method, each delay element consists of one-bit adder with the sum-out connected to the data input of a flip-flop. The two operands of the first adder are '1' and the hit signal. For the other adders, the two operands are '0' and '1'. The carry-out output of each adder is connected to the carry-in input of the following adder. Thus, before the arrival of the hit signal, i.e. the hit signal is '0', the sum-out of all the adders are set to '1'. When the hit signal arise, i.e. is set to '1', the sum-out of the first adder will be '0' and its carry-out will be set to '1'. Since the two operand of the next adder are '0' and '1', the sum-out of this adder will be '0' and its carry-out will be set to '1', etc. Hence, the hit signal will propagate through the adders' carry chain resulting zeros at the sum-out of the propagated adders, while the sumout of the un-propagated adders rest ones [5][7]. The adders' sum-out outputs are sampled by the corresponding flip-flops at the following rising edge of the system clock. Theoretically, the values of these flip-flops should be a thermometer code (unary code). But due to the clock skew and the variations in the carry delays and the routing delays from the adders to their corresponding flip-flops, there will be a bubble effect [START_REF] Xia | A 6.6 ps RMS resolution time-to-digital converter using interleaved sampling method in a 28 nm FPGA[END_REF]. An approach to overcome this problem is covered in the next section.

III. TDL TDC IMPLEMENTATION ON CYCLONE V FPGA

A. ALM structure of Cyclone V FPGAs

The logical structure of the Cyclone V is based on the logic array blocks (LABs). The LAB in Cyclone V is composed of ten adaptive logic modules (ALMs) that could be configured to perform logic, arithmetic and register functions. As presented in Fig. 2, each ALM contains four 4-input look-up tables (LUTs), two dedicated adders and four registers. The output of each adder could be sampled by two registers while the carry connections between the dedicated adders provide a fast carry chain along the LABs existing on the same column [9]. The studied TDL-based TDC makes use of the carry chain to enhance the resolution. The elementary cell of the delay line is the adder. That could be achieved by developing an M-bit adder VHDL code, proposed in [START_REF] Dadouche | Design methodology of TDC on low cost FPGA targets[END_REF], where M is the number of the delay cells. The latter is determined by the system clock period and the delay of the one-bit adder. The sum-out outputs of the cascaded adders consisting the delay line are sampled by the registers line (RL) at the rising edge of the system clock.

Two practical issues must be considered. Firstly, since the two operands of each adder are constants '1' and '0', the Quartus software would simplify the design and the delay line would not be placed when the design is compiled. To prevent that, the TDL block should be set as a design partition. Secondly, Quartus might place some adders and their corresponding registers in different ALMs, which causes a large variation in the delays from the sum-out of the adders to the data input of the corresponding flip-flops and the effect of the clock skew will be aggravated. To avoid this misplacing, both the adder block and the register block should be instantiated in an upper-level block and a logic lock region should be created for this block using the chip planner tool. The width of this region should be 1 and its height (the number of LABs consisting this region) should be calculated according to the required number of the delay elements M. Since there are 20 adders in each LAB, the region's height equals to M divided by 20. Another interest of creating a logic lock region for the TDL is to control its physical placement by specifying its origin. The effect of the physical placement on Fig. 2. ALM in arithmetic mode for Cyclone V devices [9] the TDC performance could be experimentally observed in order to choose the origin that provides the best performance.

B. Carry chain timing analyzing in Cyclone V FPGA

The timing analysis results of the carry chain path, obtained by using the Timing Analyzer tool, indicate that the propagation delay is not identical for all the adders. Furthermore, the following points were observed:

Firstly, the path report of a 1024-cell TDL shows that the propagation time from the first adder to the last one is 31.57 ns, whereas the experimental results indicates that the true propagation time is 6.4 ns. Experimental tests were performed on TDLs with different number of cells show that the timing results obtained by the Timing Analyzer should be divided by 5 to give the approximate true timing values. This fact goes along with the results provided by [START_REF] Xia | A 6.6 ps RMS resolution time-to-digital converter using interleaved sampling method in a 28 nm FPGA[END_REF]; stating that the propagation time of 691 adders is about 4 ns. This discrepancy between the experimental and the timing analyzer results is probably due to the timing analyzer model that failed to asses a good timing with such an optimized routing topology (with very short propagation delays of some ps). Indeed, the logic cell in a standard compilation are generally spread along the LAB resulting in a longer mean delay.

The second point is that the propagation delay is not identical for all the adders but differs according to the order of the adder in the ALM and the position of the ALM in the LAB. The timing analyzing results presented in Fig. 3 show that:

• The first adder existing in the first and the sixth ALM of a LAB has an approximate propagation delay of 100 ps. Dividing this value by 5 gives a real propagation delay of 20 ps for these two adders. • For the other ALMs, the first adder has an approximate propagation delay of 50 ps. Therefore, the real delay (extracted from experiments) of these adders is 10 ps. • The propagation delay of the second adder in all the ALMs is 0 ps. • The delay from the sum-out of the first adder in an ALM to the data input of its corresponding flip-flop is 425 ps (real value is 5 times lower i.e. 85 ps), whereas for the second adder, this delay is 154 ps (real value is 30.8 ps).

These differences in the propagation time of the adders as well as from the adders to their corresponding flip-flops, in addition to the clock skew, cause grate differences in the bins' widths, which leads to a nonlinear TDL. Furthermore, the TDL will have many zero-width bins and large bubbles will appear in the thermometer code at the output of the TDL. Nevertheless, the propagation time differences could be useful to improve the TDL resolution but complicated bin realignment operations should be applied to overcome the bubble effect [START_REF] Xia | A 6.6 ps RMS resolution time-to-digital converter using interleaved sampling method in a 28 nm FPGA[END_REF].

C. TDC implementation

A temporal resolution of 20 ps is generally sufficient for most TCSPC based applications that involve a Single Photon Avalanche Diode (SPAD) device as a light detector since the jitter of the SPADs is of some tens of picoseconds [START_REF] Ceccarelli | 83ps timing jitter with a Red-Enhanced SPAD and a fully integrated front end circuit[END_REF]. Hence, in order to reduce the nonlinearity and the bubble effect, and according to the timing characteristics stated above, the TDL could be designed by combining several adders in one TDL cell in a way that all the cells has the same targeted propagation times of almost 20 ps. As illustrated in Fig. 4, the following approach fulfills that:

• Each LAB is composed of an upper half and a lower half.

These two halves are identical in terms of propagation time properties; each half consists of five ALMs, and the propagation time of the ALM, comprising two adders, is (20 + 0 = 20 ps) for the first ALM and (10 + 0 =10 ps) for the next four ones.

• Three TDL cells with identical propagation time of (20 ps) could be implemented in each half LAB. The first cell is composed of the two adders of the first ALM. The second one is composed of the four adders existing in the second and the third ALMs. The third cell is composed of the four adders of the fourth and the fifth ALMs. This approach is applied for all the half LABs consisting the TDL. Fig. 5 presents the behavioral VHDL code developed to link the adders' line to the registers line for a 256-cell TDL TDC in order to apply this approach.

It is important to note that due to the high speed of the carry chain in the Cyclone V FPGA family, a very high frequency is required to clock a short TDL line if the cell is composed of only one adder (667 MHz for 256-cell TDL). Otherwise, for a reasonable frequency, a large number of oneadder cells is required to cover all the clock period (855 cells for 200 MHz) which increases the encoding time and thus the TDC's dead time. Hence, our design provides a good compromise (256 cells for 200 MHz) since the resolution is not the most critical factor for the intended application.

IV. DESIGN METHODOLOGY

A. Global system structure

Fig. 6 demonstrates the global structure of the studied TDC system. In this system, as most of the TCSPC-based systems, the light emission is piloted by a synchronous signal, called the START signal, generated by the system (the START_Generator block). Therefore, TFine1, mentioned in the timing diagram illustrated in Fig. 1, is equal to 0 and the TDC is reduced to a coarse block and one fine block [START_REF] Dadouche | Design methodology of TDC on low cost FPGA targets[END_REF], and the equation ( 1) is simplified to: The START signal, which is a clock-width pulse, clears the COARSE counter. After that, this counter starts incrementing at the system clock. When a photon is detected by the SPAD, the STOP-Generator block sets the STOP signal to '1'. At the rising edge of this signal, the value of the COARSE counter is latched giving the coarse interval (Tcoarse). The STOP signal also propagates through the TDL generating a thermometric code. This code is sampled by the registers line at the first clock rising edge to be encoded by the encoder block to a binary value (Tfine2). The WRITE_Data_Control block is a finite-state machine that controls the coarse-fine data writing to a (FIFO) memory. When a photon is detected and the STOP signal is high, this state machine waits for a certain number of clock periods before writing the data to the FIFO in order to give the encoder the required time to finish the encoding process. After writing the data, this block sends a clock-width pulse to the STOP_Generator to reset the STOP signal to '0'.

Tmeasured = N.Tclk -Tfine2 (2)

B. Timing considerations

In such systems, the routing topology has an important effect on the system behavior. As mentioned before, the TDL output should be sampled by the RL at the first clock rising edge that follows the STOP signal edge. In order to achieve that, the STOP signal is sampled at the clock rising edge by a flip-flop, and the output signal of this flip-flop is inverted and used to enable the RL flip-flops. Thus, the ENABLE signal is set to '0' after the first clock edge and disables the RL. The ENABLE signal should arrive to the RL during the clock period after the arrival of the STOP signal to the TDL. If the ENABLE signal arrives before the STOP signal, the TDL output will be sampled at the clock edge that comes before the propagation of the STOP signal resulting the code 255 at the TDL output and the bin 255 will have a large width. If the ENABLE signal is delayed more than a clock period, the TDL output will be sampled at the second following clock edge after the STOP signal. i.e. the STOP signal propagates all through the TDL resulting the code 0 at the TDL output and the bin number 0 will have a large width. In fact, mastering the arrival timing of the ENABLE signal is the most critical issue in our design. It was experimentally noticed that the ENABLE signal, with no constraints taken, arrives before the STOP signal. Furthermore, its arrival timing differs with every new compilation or when the TDL placement is changed. Therefore, to overcome this issue, the ENABLE signal is generated by sampling the STOP signal just before the first TDL cell. To achieve that, the flip-flop that generates the ENABLE signal is placed at the beginning of the RL and sample the sum-out of the first adder, which is reset to '0' once the STOP signal edge arrives. In this case, there will be no delay between the two signals. The timing analyzer results show that the ENABLE signal arrives to all the RL flip-flops in less than one clock period. Furthermore, this approach guarantees the same behavior when the design is recompiled. Nevertheless, it is recommended to perform an incremental compilation.

V. CONCLUSION

In this work, we have presented a design methodology for implementing a TDL TDC on Cyclone V FPGA targets. We have covered the ALM structures of this FPGA family and proposed an approach to overcome the delay differences along the carry chain, taking in consideration the different routing constraints that should be applied to force the desired placement of the TDL. We have also provided a simple approach that ensures the same functionality when the design is recompiled. The studied design was implemented and tested using a Cyclone V FPGA target, the achieved resolution is 19 ps which is sufficient for many TOF-based applications.

Fig. 1 .

 1 Fig. 1. Time to digital converter's timing diagram.

Fig. 3 .

 3 Fig. 3. Timing Analyzer: Carry chain path report .

Fig. 4 .

 4 Fig. 4. Chip planner: Delay cells and registers in a logic array block.

Fig. 5 .

 5 Fig. 5. Linking adders to registers line VHDL model.

Fig. 6 .

 6 Fig. 6. Design's global structure .

ACKNOWLEDGMENT

The authors would like to thank the French National Research Agency (ANR-15-CE11-0006) and the European Interreg WPS Project for funding this study.


Wassim Khaddour

Wilfried Uhring
email: wilfried.uhring@unistra.fr

Foudil Dadouche

Vincent Frick

Morgan Madec

Time-Resolved Fluorescence Measurement System for real-Time high-Throughput Microfluidic Droplet Sorting

Keywords: FPGA, Real-time, Time correlated Single Photon Counting, SPAD, Microfluidic, Droplet sorting, Screening

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

High-Throughput Screening of biomolecule in the field of pharmacology and biochemistry is a technic used to identify within a chemical library the molecules that have new properties or are biologically active. A microfluidic circuit coupled with fluorescence detection is an excellent tool for the application of High Throughput Screening (HTS) [START_REF] Miller | High-resolution dose-response screening using droplet-based microfluidics[END_REF] and Fluorescence Activated Cell Sorting (FACS) [START_REF] Baret | Fluorescence-activated droplet sorting (FADS): efficient microfluidic cell sorting based on enzymatic activity[END_REF]. Indeed, the fluorescence Quantum Yield (QY) of fluorescently labeled molecules can be significantly reduced upon interaction with other molecules. In Förster Resonance Energy Transfer (FRET)-based assays, the fluorescence yield of the excited "donor" molecule is quenched if its excitation energy is efficiently transferred to a FRET "acceptor" in close proximity. As a consequence, the measure of the fluorescence QY is the indication of the biomolecular interaction. Generally, the fluorescence QY is detected by measuring the fluorescence intensity, but this last can be affected by other parameters such as fluorophore concentration or excitation light intensity. On the contrary, Fluorescence Lifetime (FLT) detection is potentially much more accurate to assay biomolecular interactions because the latter is an intrinsic determination of the fluorescence QY independent of the other parameters. This is the reason why Fluorescence Lifetime Imaging Microscopy (FLIM) has been intensively developed for the analysis of cell biology [START_REF] Bastiaens | Fluorescence lifetime imaging microscopy: spatial resolution of biochemical processes in the cell[END_REF]. For the same reason, FLT detection is expected to enhance the reliability of biomolecular interaction assays [START_REF] Gakamsky | Use of fluorescence lifetime technology to provide efficient protection from false hits in screening applications[END_REF][5]6]. The current commercial plate readers for HTS are able to read plate up to 1536 microwells plate within a few minutes in the FLT mode [START_REF] Hansson | A Comparative Study of Fluorescence Assays in Screening for BRD4[END_REF], whereas a rate of a few thousands of droplets per second can be achieved with a microfluidic channel. Thus, the use of microfluidic can potentially accelerate the drugs discovery by a factor of 100. Nevertheless, merging both the FLT-based assays and the microfluidic is a challenging task as a droplet rate of 1 kHz requires operating both the acquisition and the lifetime assessment in less than 1 ms. Moreover, the thin thickness of about 50 µm of the droplets in the microchannel leads to a low absorption and thus a weak fluorescence signal. It has been demonstrated that the Time-Correlated Single Photon Counting (TCSPC) approach can be used to measure the FLT in microfluidic even in high throughput conditions [START_REF] Léonard | High-throughput time-correlated single photon counting[END_REF]. Nevertheless, a droplet sorting rate of only a few Hertz has been achieved up to now [START_REF] Hasan | Fluorescence lifetime-activated droplet sorting in microfluidic chip systems[END_REF].

In this paper, we present a fluorescence life-time measurement system for real-time high-throughput microfluidic droplet sorting designs. This system is implemented using a Cyclone V FPGA SoC platform. It is able to operate at high rates of more than 1 thousand droplets per second. In section II the global architecture of this system is described. Section III covers the data management and processing. The experimental results are presented in section IV. Finally, section V presents a brief conclusion.

II. SYSTEM ARCHITECTURE

The proposed system is based on the TCSPC technique. It drives a laser diode that emits a short pulse of light to excite the fluorescent substances in the passing droplet. The photons emitted by the fluorescent droplet are detected by a Single Photon Avalanche Diode (SPAD). A Time to Digital Converter (TDC) measures the arrival times of these photons. The TDC's data are then transferred to the HPS part at a high rate to be processed by a bare-metal-developed C program. This program creates a histogram that represents the fluorescence decay and measures the droplet's FLT according to which this droplet will be sorted. This section describes the architecture of the system illustrated in Fig. 1.

A. Laser diode and driving electronics

TCSPC technique requires short light pulses generated at a high repetition rate. The circuit described in [START_REF] Uhring | A low-cost high-repetition-rate picosecond laser diode pulse generator[END_REF] is used to drive a laser diode with an emitting wavelength selected according to the application, i.e. the measured fluorophore. This circuit is triggered by the system itself at an adjustable rate, ranging from a few Hertz up to 100 MHz, and allows to emit a pulse of light of typically less than 100 ps FWHM. The emitted optical mean power is about 1 mW at 80 MHz.

B. Photon detector

A commercial SPAD id100-MMF100 is used as a photon detector in this system. Its dead time is 50 ns which currently limits the maximal detectable photon rate. The use of an array of SPAD is a key to overcome this limitation [START_REF] Malass | Parallelized integrated time correlated photon counting system for high photon counting rate applications[END_REF] and represents a way to further improve the system performance.

C. FPGA Part

Time to Digital Converter (TDC):

A tapped delay line-based TDC is implemented in the FPGA part by following the methodology described in [START_REF] Dadouche | New design-methodology of high-performance TDC on a low cost FPGA targets[END_REF]. The TDC operating clock frequency is set to 200 MHz, which leads to 5 ns coarse resolution, and the tapped delay line allows to achieve a fine resolution of 19 ps.

Frequency Divider block:

The laser diode is driven by the system itself by means of the START signal. The frequency divider block generates this signal at a configurable frequency that defines the time window of the histogram. Another role of this block is to generate the PACKET signal. This signal allows the temporal packetizing of the arrival times data when transferring these data to the HPS part. The frequency of the PACKET signal is user-configurable and defines the packet duration that should be experimentally set according to the droplet size and the flow rate; a droplet should pass the detection spot during a large number of packet durations. This packetizing notion is used for the real-time droplet detection as explained later.

Write_Data_Control block:

This block is a finite-state machine that controls the writing of the arrival times data into a 32-bit FIFO before transferring these data to the SDRAM of the HPS part. When a new photon is detected and its arrival time is measured by the TDC block, the Write_Data_Control block writes the photon arrival time data to the FIFO. The TDC output is written to the 16 least significant bits of the written word as shown in Fig. 2-A. Also, when a new PACKET pulse arrives, this block writes to the FIFO a special code "0xFFFF", called the packet stamp, as illustrated in Fig. 2-B. This approach allows to packetize the arrival times data of the photons detected during the same packet duration between two packet stamps. Moreover, with every new packet pulse, an 8-bit counter is incremented, indicating the packet number that represents the absolute time of this packet. Another counter is reset by the PACKET signal and incremented with every new detected photon to count the number of photons in the packet, i.e. this counter measures the light intensity detected during the packet. These two counters are used for the data preprocessing, more precisely, for the droplet detection illustrated in Fig. 3, as follows:

• The light intensity detected during a droplet packet is much greater than that detected during a background one. Therefore, two intensity thresholds, called the droplet and the background thresholds, are experimentally defined to classify the packets. The packet is considered as a droplet packet if its photons number is greater than the droplet threshold, whereas it considered as a background packet if its photons number is less than the background threshold. Finally, a packet is unclassified if its photons number is between the two thresholds. • In order to discriminate the droplets data from the background, two one-bit flags, called the droplet and the background flags, are added to the packet word. In order to ignore the false droplet packets presenting in the background (e.g. packet X1 and X2 in Fig. 3) and the false background packets presenting in a droplet (e.g. packet Y in Fig. 3), the start of a droplet is determined by at least three successive droplet packets, i.e. the droplet flag is set to one after three successive droplet packets. This flag is reset to zero at the start of the next background, which is considered at the first background packet that comes after at least two successive unclassified packets (packet B in Fig. 3). Starting from this packet until the start of the next droplet, the background flag is set to one for every packet below the background threshold. Conversely, this flag is reset to zero for every packet above the background threshold. • A 6-bit counter is incremented at the start of every detected droplet to give this droplet a unique identifying number. • The classification flags as well as the droplet number and the packet number are added to the packet stamp to be written to the FIFO as a 32-bit packet word as shown in Fig. 2-B.

Hardware design using Intel FPGA tools

In such real-time systems, the main limitation is the data transmission and processing rates. Transferring the data to a powerful workstation is time-consuming, and the non-real time operating systems do not allow handling the high data rate of the continuous droplet stream [START_REF] Hasan | Fluorescence lifetime-activated droplet sorting in microfluidic chip systems[END_REF] [START_REF] Lieske | Embedded fluorescence lifetime determination for high throughput real-time droplet sorting with microfluidics[END_REF]. In order to overcome this limitation, a system on chip (SoC) FPGA platform is used to implement our system, that allows to process the data locally on the embedded ARM processor. For the data transmission, we designed a hardware structure, using the Intel FPGA Platform Designer integration tool, this design takes advantage of the FPGA-to-SDRAM (F2S) interface that allows the FPGA to access the SDRAM of the HPS at high transfer rates [START_REF]Modular Scatter-Gather DMA core[END_REF]. Moreover, this design also integrates a modular scatter-gather direct memory access (mSGDMA) intellectual property (IP) module that allows the data transfer by applying preloaded instructions, called descriptors, without the direct interfering of the processor in the data movement operations. Beside the HPS and the mSGDMA IPs, the hardware design involves parallel I/O (PIO) IPs that allow the HPS system to configure the FPGA system's parameters, i.e. the START and PACKET signals' frequencies, and the droplet and the background thresholds.

D. HPS software development

A bare-metal C program is developed to be run on the HPS ARM processor and perform all the data management; it configures the mSGDMA to perform the data transmission from the FPGA fabric to the HPS SDRAM, and performs the various data processing stages to calculate the droplets' FLTs according to which a sorting control signal is generated. The next section discusses the data management algorithm.

III. DATA MANAGEMENT

A. Data transfer

As mentioned before, the processor loads the mSGDMA by specific descriptors that indicate, among others, the data size and the memory addresses into which the data will be transferred. While the mSGDMA is performing the data transfer, the processor is free to perform other tasks (i.e. the data processing). This feature allows the real-time data processing procedure, as follows:

• Two memory blocks of the same size as the descriptors are allocated. • While the data is being transferred to one of these blocks, the processor processes the data already transferred to the other block as illustrated in Fig. 4.

As explained above, the mSGDMA transfers the data in predefined constant size and not in continuous stream, i.e. it should wait for a certain number of photon and packet words to be written before performing a data transmission, whereas the number of photons in a droplet (the droplet data size) is variable and unknown. In order to perform the data transfer at a rate of a droplet per transmission operation whatever the number of photons in the droplet is, the Write_data_control block writes the arrival times data into the FIFO at constant rate (the system's maximum detectable photon rate) by writing the word "0x00000000" when there is no TDC's data to be written (these zero words will be ignored when processing the data), and setting the descriptor data size to (data writing rate / droplet flow rate). Thus, when the mSGDMA is transferring a droplet's data, the processor processes the previous droplet's data, which enables the real time operation.

B. Data processing

The first step of the data processing is to create the droplets and the background histograms. Since the arrival times data are packetized with packet words tagged by the droplet and the background flags, a raw histogram for each droplet is directly created from the droplet packets data (the arrival times data comprised between Packet A and B in Fig. 3), and the background histogram is updated by adding the background packets (the arrival times data between Packet B and C). Then, a calibration routine is called to compensate the non-linearity of the TDL TDC. Thereafter, in order to eliminate the background fluorescence signal detected during the passage of a droplet, the background's mean histogram is subtracted from the droplet histogram to generate the droplet's final histogram. The latter represents the pure fluorescence decay of the droplet and is used to estimate its FLT by applying a maximum Likelihood-based algorithm [START_REF] Lieske | Embedded fluorescence lifetime determination for high throughput real-time droplet sorting with microfluidics[END_REF]. According to the measured FLT, an action could be taken to sort the droplet. The droplets' information, including the droplet's number and the packet number of its first packet as well as an action-bit, which is set or cleared according to the measured FLT and the sorting criteria, are written to a FIFO in the FPGA part. According to these information, the Sorting control block generates the CONTROL signal that drives the actuator to perform the droplet sorting.

IV. EXPERIMENTAL RESULTS

A. Experimental setup

This system was implemented using a low cost SoC FPGA platform, namely the "Altera DE-10 Nano" kit. As a pulsed light source, a 450 nm, 50 mW CW laser diode is driven by a fast pulse generator described in [START_REF] Uhring | A low-cost high-repetition-rate picosecond laser diode pulse generator[END_REF]. The laser diode is mounted on an optical lens system to collimate the beam. For the photon detection, an ID100 commercial SPAD from IDQ is used and coupled with an optical 480 nm high pass filter. In order to test our system, the droplet-background sequence was emulated using some white paper segments separated by nonfluorescent dark parts. The set rotates at a constant speed thanks to a rotating disc acting as an optical chopper. The droplets' sizes are defined by the segment width, while the flow rate is determined by the chopper's rotation rate; the droplet flow rate is calculated by multiplying the number of paper segments by the motor rotation rate. Fig. 5 presents the system set-up. Most of the system's parts cost less than 200$ each except for the commercial SPAD.

B. Experimental results discussion

In order to evaluate the system performance, extensive experimental tests were performed at different flow rates and with droplets of different sizes. The obtained results showed that the system was able to detect all the emulated droplets without reporting any missed droplet. Since all the droplets have the same FLT because they are made of the same material, the measured FLT values has a standard deviation close to the theoretical quantum photon shot noise limit, even for a low number of photons in droplet. Table I shows the experimental results for three emulated flow rates, whereas Fig. 6 shows the droplet histogram at these flow rates. The obtained results show that even at a high flow rate of 1000 droplet per second, the signal to noise ratio remains above 40 allowing efficient screening as explained in [START_REF] Léonard | High-throughput time-correlated single photon counting[END_REF]. V. CONCLUSION

In this work, we have presented a FLT measurement system for real-time high-throughput microfluidic droplet sorting. We have described the global structure of this system and introduced the packetizing notion, which provides a robust approach for droplet detection. We have tested this system on an emulated droplet flow at different rates (up to one thousand droplets per second) and with different droplet sizes. The system detected the droplets and measured their FLT. The standard deviation of the measured FLT values is only 30% above the theoretical quantum photon shot noise limit. This system paves the way to FLT droplet sorting in high-throughput condition with microfluidic.

Fig. 2 .Fig. 1 .

 21 Fig. 2. Photon and Packet words. A-Photon arrival time data word, B-Packet word.

Fig. 3 .

 3 Fig. 3. Droplet detection: 2 thresholds are set to detect the start and the end of the droplet with a robust algorithm insensitive to noise spikes.

Fig. 5 .

 5 Fig. 5. Experimental set-up. 1-SoC FPGA platform with power supply, 2-Pulsed laser driver, 3-Laser diode with collimating lens, 4-480 nm high pass optical filter, 5-SPAD, 6-Optical Chooper used to simulate the droplets flow, 7-Chooper controlor.

Fig. 4 .

 4 Fig. 4. Data transfert flux diagram.

TABLE I . EXPERIMENTAL RESULTS Flow rate (droplet/ sec) Average number of photons in droplet Fluorescence lifetime (FLT) values Average (ps) Experimental Theoretical Standard deviation (ps) SNR Photon shot noise limit (ps) SNR

 I

	5	6,39E+05	973	1.6	617.4	1.2	804.1
	200	1,63E+04	976	10	97.1	6.8	144.6
	1000	3,20E+03	970	22	43.5	17	56.7

ACKNOWLEDGMENT

The authors would like to thank the French National Research Agency (ANR-15-CE11-0006) and the European Interreg WPS Project for funding this study.


Suhail Mohd
email: mohd-suhail.rizvi@univ-grenoble-alpes.fr

Rizvi

Abdessamad Nait-Ouhra

Alexander Farutin

Philippe Peyla

Salima Rafai

Chaouqi Misbah
email: chaouqi.misbah@univ-grenoble-alpes.fr

Rheological signature of microswimmer phase-locking under flow

Understanding macroscopic behaviors of suspensions of active particles is a challenging issue in the study of living fluids. We investigate the response of shape-change driven microswimmers to an external flow. In particular we estimate the viscosity of the dilute suspension of active microswimmers for the flow timescale close to that of the swimmer activity. We find that the shape activity leads either to an increase or decrease of the effective viscosity, depending on the shear rate. These opposite behaviors occurring for the same microswimmer originate from a phase-locking phenomenon between the swimmer shape oscillations and the applied flow. A simplified analytical model remarkably reproduces the results of full simulations, and offers a promising framework for the derivation of macroscopic constitutive laws for active matter.

I. INTRODUCTION

Because of the wide variety of the environments and ecosystems where micro-organisms live, their evolution has resulted in various ways of moving (gliding, twitching, swimming) towards their favorite chemical or food sources [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF][START_REF] Elgeti | Physics of microswimmers-single particle motion and collective behavior: a review[END_REF]. Examples of propulsion machinery include the cyclic beating of two flagella in a breast-stroke manner by C. reinhardtii [START_REF] Elgeti | Physics of microswimmers-single particle motion and collective behavior: a review[END_REF], rotation of helical flagellar bundle by E. coli [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF] or the repeated body-deformation of the Eutreptiella gymnastica [START_REF] Throndsen | Flagellates of Norwegian coastal waters[END_REF]. Similarly, cells of the immune system develop ample shape changes to fight pathogens in complex environments, such as the extracellular matrix within tissues. There is now increasing evidence that cells of the immune system, among others, are able to swim in a fluid [START_REF] Barry | Dictyostelium amoebae and neutrophils can swim[END_REF][START_REF] Bergert | Force transmission during adhesion-independent migration[END_REF]. The micron size of cells and their low velocity ensure that inertia is negligible. Despite the diversity of micro-organisms, they all share a common feature, namely the repetitive and non-reciprocal temporal variation of forces on the fluid.

A popular swimmer model is the constant force-dipole which often discards the details of the active force generation machinery and retains the average behavior over one cycle of time-dependent active forces [START_REF] Saintillan | The dilute rheology of swimming suspensions: A simple kinetic model[END_REF][START_REF] Drescher | Direct measurement of the flow field around swimming microorganisms[END_REF]. This swimmer can be categorized into two types-"puller" (eg. Chlamydomonas reinhardtii [START_REF] Drescher | Direct measurement of the flow field around swimming microorganisms[END_REF][START_REF] Rafai | Effective viscosity of microswimmer suspensions[END_REF]) and "pusher" (eg. E. coli [START_REF] Drescher | Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering[END_REF]) depending on the sign of the dipole. A suspension of microswimmers displays unusual macroscopic rheology [START_REF] Rafai | Effective viscosity of microswimmer suspensions[END_REF][START_REF] Hatwalne | Rheology of active-particle suspensions[END_REF][START_REF] Saintillan | Extensional rheology of active suspensions[END_REF][START_REF] Takatori | Superfluid Behavior of Active Suspensions from Diffusive Stretching[END_REF][START_REF] Gachelin | Non-Newtonian viscosity of Escherichia coli suspensions[END_REF][START_REF] Lopez | Turning Bacteria Suspensions into Superfluids[END_REF]: a pushers' (pullers') suspension exhibits a lower (higher) effective viscosity than their respective passive counterparts. Furthermore, experimental studies have also shown a rich rheological behavior rarely recovered in theoretical models. For instance, a shear-thinning [START_REF] Rafai | Effective viscosity of microswimmer suspensions[END_REF] (shear-thickening [START_REF] Gachelin | Non-Newtonian viscosity of Escherichia coli suspensions[END_REF][START_REF] Lopez | Turning Bacteria Suspensions into Superfluids[END_REF]) nature of the pullers' (pushers') suspensions has been reported. So far, the main key ingredient of the models explaining the rheology of microswimmers (increase and decrease of viscosity due to activity), is rotational diffusion [START_REF] Saintillan | Rheology of Active Fluids[END_REF]. However, in general, microswimmers exhibit shape deformations (flagella motion, body shape changes, etc.) and their impact needs to be explored. It has been shown recently that the cyclic motion of flagella makes the cell follow very specific Jefferey's orbits that can affect rheology [START_REF] Jibuti | Effective viscosity of a suspension of flagellar-beating microswimmers: Three-dimensional modeling[END_REF].

In order to decipher the rheological behavior of swimmer suspensions, it seems relevant to consider the coupling between the flow and the time-dependent active forces. More generally, in most cases, shape changes of microswimmers are intrinsically related to propulsion mechanisms (flagellar beating or shape deformations). We show that shape change activity not only contributes to the viscosity of the suspension without evoking noise, but also uncovers new subtle effects. In particular, the active contribution to the viscosity is found to switch between positive and negative values depending on the magnitude of the applied shear flow for a given type (puller or pusher) of microswimmer. This results from a nonlinear coupling between the swimmer and the applied flow.

A widely adopted model for swimmers consists of a thin rod-like rigid shape with a given stresslet (force dipole), which is aligned with the rod and is either inward (puller) or outward (pusher). A deterministic rod under linear shear flow is known to exhibit tumbling (Jeffery orbit) and spends the same time in the upper and lower halves of the shear plane. As a consequence, no resulting active contribution to the viscosity is expected on average. Noise breaks this symmetry (the rod, then, spends more time in the upper half plane), giving rise to a net contribution of the active stress to the rheology [START_REF] Saintillan | The dilute rheology of swimming suspensions: A simple kinetic model[END_REF][START_REF] Takatori | Superfluid Behavior of Active Suspensions from Diffusive Stretching[END_REF]. Hence, in all available models, noise (due, for instance, to rotational diffusivity or tumbling of the microswimmers) is a crucial ingredient to recover a net effect of average active stresses on viscosity. Real swimmers undergo repeated shape changes, that lead to a non-trivial contribution to the active stress. Here we will address the effect of shape activity on the rheology of microswimmer suspensions in two cases: amoeboid swimmer and bead-spring microswimmer.

II. RHEOLOGY OF MICROSWIMMER SUSPENSIONS

A. Amoeboid microswimmer

We begin our study by considering an amoeboid microswimmer in 2d-space (Fig. 1) characterized by a fluid-filled impermeable membrane of perimeter P 0 and enclosed area A 0 immersed in another fluid [START_REF] Wu | Amoeboid motion in confined geometry[END_REF][START_REF] Wu | Amoeboid swimming in a channel[END_REF]. The swimmer has reduced area Γ = 4πA 0 /P 2 0 (measuring deviation from a circle for which Γ = 1) which allows the microswimmer deformation. η and λη denote the respective viscosities of external and internal fluids. We consider the flow, in and outside of the swimmer to be in Stokes regime, following

∇ • u = 0, η i ∇ 2 u (r, t) = ∇P (r, t) (1)
where u and P are the fluid velocity and pressure, respectively, and η i = η or λη depending on the location of r.

The shape of the swimmer at time t is described by a closed curve X(s, t) in x-y plane, where s is the curvilinear coordinate along the membrane. We represent the membrane position X by a Fourier series

X(s) = kmax k=-kmax X k exp (2πiks) (2)
where the complex amplitudes X k are the shape parameters of the membrane. Given the inextensibility of the membrane the application of active force results in generation of tension ξ(s, t) which gives rise to a tangential force ξ(s, t)t(s, t). The resultant passive force density on the swimmer surface is, therefore, given by

f (s, t) = -ξ(s, t)c(s, t)n(s, t) + dξ(s, t) ds t(s, t) (3)
where n and t are the normal and tangential unit vectors in the membrane surface, respectively. In addition to the passive membrane forces, there is also an active force normal to the membrane which is also written in terms of Fourier harmonics as

F a (s, t) = n=nmax n=-nmax,n =-1,0,1 f n (t)e ins . (4)
In the above, due to the fluid incompressibility the term corresponding to n = 0, which corresponds to a constant pressure jump across the membrane, plays no role. For simplicity we set f ±1 (t) = 0 (we have kept the model as simple as possible by considering only the second and third harmonics). We require at least two harmonics to induce shape deformations which are not invariant under time-reversal a necessary condition for autonomous swimming (Scallop theorem [START_REF] Purcell | Life at low reynolds number[END_REF]). We set f 2 = -f a cos(ω a t)/2 and f 3 = f a sin(ω a t)/2, where ω a is the stroke frequency and f a is the active force amplitude. By considering only the real part of the Fourier harmonics expansion, we get the following active force

F a (s, t) = -f a cos(ω a t) cos (2s) + f a sin(ω a t) cos (3s) (5)
Including higher harmonics is not decisive for the swimming phenomenon [START_REF] Farutin | Amoeboid swimming: a generic self-propulsion of cells in fluids by means of membrane deformations[END_REF]. Therefore, the total force density on the swimmer surface is given by

F (s, t) = F a (s, t) n + f (s, t) + f o (t) + f t (t)t(s, t) (6)
where f o (t) and f t (t)t(s, t) allows to fulfill the force-free and torque-free conditions on the swimmer surface. We define the non-dimensional active force amplitude f = f a / (ηω a ) as the ratio between the time scales of the swimming strokes and fluid flow. This, together with hydrodynamic interactions, results in autonomous propulsion along a straight line [START_REF] Wu | Amoeboid motion in confined geometry[END_REF][START_REF] Farutin | Amoeboid swimming: a generic self-propulsion of cells in fluids by means of membrane deformations[END_REF][START_REF] Wang | Computational analysis of amoeboid swimming at low Reynolds number[END_REF]. We make use of Green's function techniques to write the membrane velocity as an integral equation expressing non-local hydrodynamics effects [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF]. Details of the boundary integral method based numerical technique are given in the Appendix-A. The system is subjected to a linear shear flow of strength γ (Fig. 1) and the effective viscosity η ef f = σ xy / γ is measured as the ratio of average shear stress (active+passive) and applied shear rate. For a suspension with swimmer volume fraction φ, the intrinsic viscosity of the suspension [η] = (η ef f -η) /φη is obtained by the Batchelor's formula [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF] as

[η] = 1 A 0 γη -F x yds + η(λ -1) (u x n y + u y n x )ds (7)
where x and y subscripts denote the x and y-components respectively, and u is the membrane velocity. Here the integral is performed over the membrane and • refers to time average (to be specified later). We define the contribution of activity to viscosity as ∆η = [η] -[η] p where [η] p is the intrinsic viscosity of the passive suspension, obtained by setting the active force to zero.

The dependence of intrinsic viscosity on the applied shear rates is summarized in Fig. 2. We observe that the contribution of activity switches from positive to negative at some critical shear rate γc . Three successive regimes of [η] are observed. For low values of shear rates ( γ < γ-), d [η] /d γ > 0 followed by d [η] /d γ < 0 for intermediate values of shear rates. Finally for high shear rates ( γ > γ+ ), [η] again starts to increase with γ. The transitions occur at shear rates γ∓ which depend on the active force amplitude. A remarkable observation is that the normalized transition shear rates ( γ∓ / γc -1) /f do not depend on the active force amplitudes f (Fig. 2B) (see Sec.-III A for the motivation for this scaling.)

In 2D a dilute passive suspension of disks provides [η] p = 2 (equal to 5/2 in 3D) [START_REF] Brady | The einstein viscosity correction in n dimensions[END_REF][START_REF] Doyeux | Effective viscosity of two-dimensional suspensions: Confinement effects[END_REF]. For vesicles, it is about 1.8 for Γ = 0.75 (Fig. 2A,black). In the presence of activity, the intrinsic viscosity attains values close to 3 (see Fig. 2A, green), showing the importance of activity (in the same proportion as found experimentally in 3D for Chlamydomonas [START_REF] Rafai | Effective viscosity of microswimmer suspensions[END_REF]). From data shown in Fig. 2A (inset) obtained for two different deflation parameters Γ we expect this value to be significantly larger than 3 for more deflated shapes such as that of Eutreptiella gymnastica. An estimate of active forces generated by the cells shows that the dimensionless force f 1 (assuming η ≈ 1cP and ω a ≈ 1s -1 ) [START_REF] Fouchard | Three-dimensional cell body shape dictates the onset of traction force generation and growth of focal adhesions[END_REF]. Owing to the membrane incompressibility [η] saturates beyond a value of f of few tens for Γ = 0.75 (Fig. 2B). In order to gain insights into this rheological behavior, we estimated the effective viscosity of the dilute suspension of a much simpler system of active swimmer composed of beads and springs.

We consider bead-spring microswimmer where three identical beads (radius a each) are connected by linear springs (length l a, stiffness k = force/strain) in a triangular manner (Fig. 1) [START_REF] Rizvi | Three-bead steering microswimmers[END_REF]. This is a variant of the model of flagellar swimming of C. reinhardtii proposed in [START_REF] Friedrich | Flagellar synchronization independent of hydrodynamic interactions[END_REF][START_REF] Kyriacos | Antiphase synchronization in a flagellar-dominance mutant of chlamydomonas[END_REF]. In the Stokes regime, that is Re = 0, we can write the equations of motion for each bead as

dr i dt = µf i (r, t) + j =i G ij • f j (r, t) + u ∞ (8)
where µ = (6πηa) -1 is the bead mobility, u ∞ is the externally imposed flow and G ij is the Oseen tensor

G ij (r ij ) = 1 8πη I r ij + r ij ⊗ r ij r 3 ij (9)
describing the hydrodynamic interaction between the ith and jth beads, and f i is the total force acting on ith bead. By virtue of the hydrodynamic interactions among the beads, the active forces f ij = f a sin(ω a t + α ij ) between ith and jth beads along the connecting spring lead to autonomous translation and/or rotation for appropriate values of α ij [START_REF] Rizvi | Three-bead steering microswimmers[END_REF]. In the absence of any external flow, for α 12 = α 13 = 0 and α 23 = α = 0 the swimmer performs translational motion. In general, the nonlinear nature of Eq. ( 8) makes it almost impossible to solve analytically. However, for f a k, the perturbation analysis around the equilibrium shape of the triangular microswimmer makes Eq. ( 8) analytically tractable. It has to be noted that the stresslet associated with this triangular microswimmer vanishes over the course of one cycle of active force. Therefore, it cannot be classified as puller or pusher microswimmer [START_REF] Rizvi | Three-bead steering microswimmers[END_REF].

Under a linear shear flow, the contribution of the triangle to the suspension viscosity is given by the Kramers-Kirkwood formula [START_REF] Doi | The Theory of Polymer Dynamics[END_REF] [η] = -

1 v b γη i y i (t)f x i (t) (10)
where v b is the swimmer volume, y i (t) is the y-coordinate of the ith bead and f x i is the x-component of the total (active+passive) force acting on it. The resulting rheology is obtained numerically and shows the same qualitative features as in the case of amoeboid model (Fig. 3). To unravel the mechanisms that govern the observed transitions (from positive to negative active viscosity as a function of shear rate) we will, in the following, analytically quantify the coupling of the active forces and applied flow.

III. ANALYTICAL CALCULATIONS A. Triangle dynamics and phase locking

For a passive vesicle or bead-spring triangle, one can describe the respective dynamics under shear flow in terms of the particle orientation ϕ s (see Appendix-C for details). Therefore, the time-averaging in Eqs. ( 7) and [START_REF] Hatwalne | Rheology of active-particle suspensions[END_REF] for the passive particles can also be achieved by an averaging over ϕ s . On the other hand, activity affects shape and thus the swimmer orientation. The phase related to activity is ϕ a = ω a t, whereas the phase related to shear flow driven shape change is ϕ s = ω s t where, due to symmetry, the shape changes occur at a frequency twice of the tumbling frequency (ω s γ; as for the classical Jeffery orbit). This means that the shape configuration depends on both phases. For γ ω a and γ ω a (weak coupling), ϕ s can be considered independent of the activity and the time-averaging in Eqs. [START_REF] Drescher | Direct measurement of the flow field around swimming microorganisms[END_REF] and [START_REF] Hatwalne | Rheology of active-particle suspensions[END_REF] has to be performed over the larger period (4π/ γ or 2π/ω a ). For γ ≈ ω a , however, there is a strong coupling between the swimmer orientation ϕ s and its activity ϕ a . For weak deformation we have (which is an exact relation for rigid object)

dϕ s dt ≈ 1 l 2 3 i=1 r i (ϕ s , ϕ a ) × ṙi (ϕ s , ϕ a ) ( 11
)
where r i is the position vector of the ith bead. The evolution equation of a bead position is expressed as a proportionality between ṙi (ϕ s , ϕ a ) and the total force acting on the bead. The active part contains obviously ϕ a and the beads positions (via the spring force), which can be written as a function of ϕ s to leading order in deformation. Together

dδ dt ≈ -( γ -ω a ) + A γ sin(δ + δ 0 ) (12)
where δ = ϕ s -ϕ a ,

cos δ 0 = 3 -3 cos α -2 sin α/ζ 2 (1 -cos α) (9 + 4/ζ 2 ) . ( 13
)
and

A = 2 √ 2 √ 1 -cos α 9 + 4/ζ 2 f a k (14)
depends on the magnitude of the triangle deformation (f a /k) due to the active force (see Appendix-B for more general form of the Adler equation for active triangle). In the above expressions, ζ = kµ/ω a l is the ratio of the timescales associated with active force (∼ 1/ω a ) and elastic deformation (∼ l/kµ). Equation ( 12) is a well known equation in nonlinear systems and is used as a model for synchronization of oscillators [START_REF] Kuramoto | Self-entrainment of a population of coupled non-linear oscillators[END_REF][START_REF] Pikovsky | Synchronization, A universal concept in nonlinear sciences[END_REF]. It is apparent from the Eq. ( 12) that for -A ≤ ( γ/ω a -1) ≤ A there is a "locking" of the two phases as dδ/dt = 0 and we obtain γ∓ = ω a (1 ∓ A) as the boundaries of the "phase-locked" regime (see Movies in [START_REF]γ-< γ < γ+ and γ > γ+) in the reference frame co-moving the with the swimmer center of mass[END_REF]). For the amoeboid swimmer also we have the shear rates (γ ∓ ) proportional to the active force amplitude (Fig. 2B), a scaling inspired by Eq. [START_REF] Takatori | Superfluid Behavior of Active Suspensions from Diffusive Stretching[END_REF]. The phase locking phenomenon can be appreciated by analyzing the dependence of angular frequency of the shear flow driven triangle dynamics (ω s ) on the applied shear. By definition the angular frequency ω s is twice the tumbling frequency due to the symmetry of the shear flow. Fig. ( 4)A shows that ω s increases with the shear rate monotonically, except in the phase-locked regime where the synchronization of shear and activity modes results in ω s = ω a . The effect of phase locking is further reflected in the time dependence of the shear stress applied by the swimmer on fluid (Fig. 5). In the absence of activity under flow the particle has a cyclic motion with frequency ω s (with the leading deformation mode behaving as e iωst + c.c.). The mode associated with activity is e iωat + c.c.. Due to nonlinear coupling, this generates new modes having frequencies ω s ± ω a . In the phase-locked regime (ω s ≈ ω a ) the system shows two frequencies: the basic one ω a and 2ω a (Fig. 5, middle panels). Outside the phase locked regime the system shows three disparate basic frequencies (generically, not rationally related to each other) ω a and ω s ± ω a , as seen in Fig. 5 (top and bottom panels).

B. Effect of phase-locking on effective viscosity

In the phase locked regime by using Eq. ( 10) we get ∆η

[η] p ≈ [3ζ (cos α -sin α -1) + 2 (cos α + sin α -1)] (9ζ 2 + 4) (1 -cos α) ω a γ 2 γ γc -1 (15)
where γc = ω a 1 + 3ζ 2 2(1 -cos α) 9ζ 2 + 4 cos α + sin α -1 3ζ (cos α -sin α -1) + 2 (cos α + sin α -1)

f a k . (16)
This shows that depending on α, γc can increase or decrease with f a . It is remarkable that ∆η/ [η] p is independent of f a in the phase-locked regime (Fig. 3B) which is also observed for amoeboid swimmer (Fig. 2B). The effect of active force amplitude f a , however, does reflect in the maximum change in the viscosity, that is |∆η/[η] p | max which is obtained by substituting γ = γ∓ in Eq. [START_REF] Saintillan | Rheology of Active Fluids[END_REF]. Since γ∓ scale linearly with f a (at least for f a k) we have |∆η/[η] p | max ∼ f a . It needs to be noted that this linear relation between maximum viscosity change and f a holds for f a k. For large f a the deformation of the triangle with linear springs can be very large and the results based on perturbation analysis may not hold.

Further, note that higher order phase-locking may also take place and Eq. ( 12) can be extended for δ mn = mϕ s -nϕ a where m = 0 and n = 0 are arbitrary integers. The occurrence of higher order phase locking suggests multiple values of critical shear rates γc where transition between shear-thickening and shear-thinning behaviors can be observed (Fig. 6). A noteworthy point is that hydrodynamic interaction is not essential for the reported feature in Fig. 3A, albeit necessary for the self propulsion. The phase locking is found to be robust against several perturbations. For the microscopic swimmers noise has a significant effect on their behavior. In order to study the influence of noise on the swimmer behavior under shear flow and the phase-locking phenomenon we consider the noise to be of two origins. First, the noise intrinsic to the swimmer coming from the mechanism of its activity. Its intrinsic nature ensures it does not contribute to the total force and torque. Second, the noise experienced by the swimmer due to the fluctuations in the surrounding fluid flow. These fluctuations in fluid flow result in additional velocities to each bead of the triangular swimmer which we assume to be uncorrelated over time as well as among the beads. We consider the combined effect of noise in the form of random forces (proportional to the velocity fluctuations in the fluid) χ i (t) on each bead such that they are uncorrelated, that is χ i (t)χ j (t ) = |χ| 2 δ ij δ(tt ) where |χ| is the noise amplitude. This nature of the noise results in a nonzero contribution to the total force and torque to the swimmer, resulting in translational and rotational diffusion of the swimmer. We performed numerical simulations of the bead-spring swimmer with noise under shear flow using Euler-Maruyama method for the stochastic ordinary differential equations [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. As shown in Fig. 3A the noise does not affect the qualitative features of the suspension viscosity as long as |χ| f a . The consideration of noise reflects in the form of an additional diffusive term in the Adler equation (Eq. ( 27)) resulting in the diffusion of the phase-difference δ on a tilted washboard potential with occasional phase slips [START_REF] Pikovsky | Synchronization, A universal concept in nonlinear sciences[END_REF] as opposed to its deterministic behavior.

IV. MECHANISMS LEADING TO THE VISCOSITY BEHAVIOR

The microswimmers are usually categorized as puller or pusher depending on far-field velocity field generated by the swimmer. The leading order contribution in the far-field velocity field due to swimmer can be written in terms of the first moment of the force distribution, that is Σ ij = 3 i=1 f i ⊗ r j for the bead-spring swimmer where r i and f i are the position and the total force on the ith bead. In the absence of any external flow (such as shear flow in this paper), the symmetric nature of the swimmer results in the average velocity field over one active force cycle to be proportional to the average stresslet Σ s = Σ xx -Σ yy where x is the direction of the swimmer propulsion. The puller and pusher nature of the swimmer is characterized by the sign of the stresslet Σ s (Σ s < 0 for puller and Σ s > 0 for pusher). The bead-spring triangular swimmer with three identical beads is known to demonstrate a puller and pusher characteristic within an active force cycle while on average its behavior is neutral (Σ s = 0) [START_REF] Rizvi | Three-bead steering microswimmers[END_REF]. In the presence of the shear flow the triangle demonstrates tumbling and we can define the stresslet in the reference frame attached to the swimmer as Σ = Σ x x -Σ y y where the swimmer propulsion direction x is co-rotating with it. We calculated the average value (over sufficiently long time to capture the effects of interaction between shear and activity modes) of stresslet Σ of the bead-spring swimmer under different shear rates. Fig. 4B shows the transition of swimmer nature from pusher (Σ < 0) to puller (Σ > 0) at shear rate close to the angular frequency of active force. This transition between the two swimmer states coincides with the transition from ∆η > 0 to ∆η < 0. This demonstrates that for shape-deformation driven swimmers the puller-pusher nature is not an intrinsic property but depends on the external conditions, such as external flow (as shown here) or confinement [START_REF] Wu | Amoeboid motion in confined geometry[END_REF].

In order to dig further into the origin of the viscosity change, we investigated the probability density of the microswimmer surface (or beads in case of bead-spring microswimmer) being at any location in the reference frame co-moving with the swimmer center of mass. This density represents the average microswimmer shape in the steady state (Fig. 7). For low γ the average shape first has an increasing cross section against the flow (meaning shear thickening), then a decrease of cross section until alignment with the flow in phase-locked regime, then an increase of cross section with γ.

V. DISCUSSION

The phase locking is very generic due to nonlinear coupling of two oscillators. Shear flow can lead to tumbling (which is the first autonomous oscillator). This is the case, for example, for Chlamydomonas which undergoes (more or less complex) tumbling; see [START_REF] Rafai | Effective viscosity of microswimmer suspensions[END_REF]. On the other hand the activity is materialized by a periodic beating of flagella (the second oscillator). The flow affects the motion of flagella, which in turns affects the tumbling dynamics, resulting in an effective coupling between the two oscillators (the activity mode and the tumbling one). We thus expect phase-locking to occur for various swimmers since tumbling and shape activity are quite common. Phase locking (at least the first order) seems to be absent for shakers (Eq. ( 27)). It needs to be pointed out that the mechanism of viscosity modulation by the swimmer activity shown here is fundamentally different from the previously known mechanisms which rely primarily on the rotational diffusion [START_REF] Saintillan | The dilute rheology of swimming suspensions: A simple kinetic model[END_REF][START_REF] Hatwalne | Rheology of active-particle suspensions[END_REF][START_REF] Saintillan | Rheology of Active Fluids[END_REF] or an intrinsic tumbling mechanism [START_REF] Takatori | Superfluid Behavior of Active Suspensions from Diffusive Stretching[END_REF][START_REF] Saintillan | Rheology of Active Fluids[END_REF] of the microswimmer. In general, however, the two mechanisms-deterministic (shape activity as reported here) and stochastic (say diffusion or run-and-tumble) can be visible under different applied shear rates for the same microswimmer. For shear rates very different from that of the angular frequency ω a of the swimming strokes it is the stochastic mechanism which has significant contribution as opposed to the shear rates close to ω a where the deterministic mechanism takes over. Equation [START_REF] Saintillan | Rheology of Active Fluids[END_REF] shows that in the phase locked region (but for γ = ω a ) the active contribution to viscosity scales as ∆η/ [η] p ∼ (ω a / γ) 2 ∼ 1. On the other hand for fixed shapes (slender body) with noise it is of order of ∆η/ [η] p ∼ t r /15t s [START_REF] Saintillan | Rheology of Active Fluids[END_REF], where t r and t s are reorientation and swimming time (over swimmer size), respectively. For several swimmers [START_REF] Rafai | Effective viscosity of microswimmer suspensions[END_REF][START_REF] Gachelin | Non-Newtonian viscosity of Escherichia coli suspensions[END_REF][START_REF] Saintillan | Rheology of Active Fluids[END_REF] the latter leads to ∆η/ [η] ∼ 1, meaning that shape activity contribution may be of the same order as that of noise.

Since most of the experimental measurements of the active suspension viscosity has been for γ ω a [START_REF] Rafai | Effective viscosity of microswimmer suspensions[END_REF][START_REF] Gachelin | Non-Newtonian viscosity of Escherichia coli suspensions[END_REF][START_REF] Lopez | Turning Bacteria Suspensions into Superfluids[END_REF] the validation of the phenomenon reported in this paper would require further experimentation, especially for γ ≈ ω a .

For flagellar microswimmers such as E. coli and C. reinhardtii ω a ≈ 100πs -1 . Thus high enough shear rates are required for the phase-locking. In contrast, amoeboid swimmers [START_REF] Throndsen | Flagellates of Norwegian coastal waters[END_REF][START_REF] Barry | Dictyostelium amoebae and neutrophils can swim[END_REF] are endowed with much longer stroke duration, approximately few seconds, requiring much lower shear rates. Finally, it will be an interesting task for future research to study multiple swimmers in order to analyze the interplay between their mutual coupling (which may trigger their synchronization) and the effect of external flow.

B. Bead-spring microswimmer

Passive triangle in shear flow

We can obtain the passive counterpart of the active bead-spring microswimmer by setting f a = 0. Once it is placed in shear flow γ it undergoes rotation and also oscillatory deformations. For γ kµ/l it is easy to see that the time-scale of triangle rotation and deformation is 4π/ γ. For small γ we can define the phase of the passive triangle ϕ s as the angle made by one the triangle beads with some pre-specified direction in the triangle plane. Therefore, all the properties of the passive triangle under shear can be written in terms of ϕ s . We can write the contribution of the passive triangle to total shear stress in the form of Kramers-Kirkwood stress [START_REF] Doi | The Theory of Polymer Dynamics[END_REF] as

σ xy KK = 3 i=1 f x i (ϕ s ) y i (ϕ s ) ( 21
)
where averaging is performed over ϕ s ∈ [0, 2π], that is all triangle configurations. For small triangle deformation, the force on any bead can be approximated as f x,y i (ϕ s ) = 3 n=1 j =i P x,y ij sin(nϕ s ) + Q x,y ij cos(nϕ s ) , where the nϕ s dependence is due to the three-fold symmetric shape of the triangle. Similarly, the small deformation assumption also gives y i (ϕ s ) = R i sin ϕ s + S i cos ϕ s . The coefficients P x,y ij , Q x,y ij , R i and S i are determined from the equations of motion of the beads (Eq. 8). Substituting the values of these coefficients in the Eq. ( 21) gives the effective viscosity of the passive bead-spring triangle suspension to be

η p = η + φ v b σ xy KK γ ≈ η + φ v b 9l 2 4µ 1 c 2 + 9 ( 22
)
where v b is the swimmer volume, φ is the swimmer volume fraction in the suspension and c = γl/kµ is the capillary number. This demonstrates a shear-thinning behavior of the passive triangle under shear flow.

Active triangle in shear flow

a. No coupling between the shear flow driven dynamics and activity: For the active triangle in shear flow if we do not consider any coupling between the shear flow and activity driven dynamics the triangle undergoes tumbling at a constant rate (due to shear flow) and deformation (due to activity as well as flow). In the absence of the coupling, the two modes-activity and shear flow are independent of each other and all dynamical variables can be written in terms of ϕ s = γt/2 and ϕ a = ω a t. Therefore, we can follow the same steps as for the passive triangle above. However, there are additional terms due to the triangle activity in the expressions for the force on the bead and but with additional dependence on activity which gives f x,y i (ϕ s ) = j =i P a ij sin(ω a t) + Q a ij cos(ω a t) sin(ϕ s ) + P a ij sin(ω a t) + Qa ij sin(ω a t) cos(ϕ s ) + passive force, where coefficients P a ij , Q a ij , P a ij and Qa ij can be obtained by solving the Eq. ( 8) for an isolated active triangle without any external shear flow. We substitute these values in Eq. ( 21) and take an average over one cycle of active force to obtain which has a monotonic dependence on γ contrary to the observation in Fig. 2. Therefore, the interaction between the two modes of triangle deformation cannot be ignored and the assumption of ϕ s ∼ t breaks down.

η a ≈ η p - φ v b f a k 2 l 2
b. With coupling between the shear flow driven dynamics and activity: In order to study the effect of coupling we first need to obtain the dependence of ϕ s on t. Under the assumption of small deformation of the triangle, we can write

dϕ s dt ≈ 1 l 2 3 i=1
x i (ϕ s , ϕ a ) dy i dt (ϕ s , ϕ a ) -y i (ϕ s , ϕ a ) dx i dt (ϕ s , ϕ a ) [START_REF] Doyeux | Effective viscosity of two-dimensional suspensions: Confinement effects[END_REF] where (x i , y i ) are the coordinates of the i th bead. In this condition we can take x i (ϕ s (t), ϕ a (t)) = R x i sin ϕ s (t) + S x i cos ϕ s (t) + R ax i sin ϕ a (t) + S ax i cos ϕ a (t) and y i (ϕ s (t), ϕ a (t)) = R y i sin ϕ s (t) + S y i cos ϕ s (t) + R ay i sin ϕ a (t) + S ay i cos ϕ a (t). Substituting these values in the following equation for the phase difference

dδ dt = dϕ s dt - dϕ a dt (26)
gives us the time dependence of the ϕ s since the phase of the activity ϕ a = ω a t is already known. In a more general setup, we can also have the activity described not by explicit time dependence but triangle configuration. We do not consider those cases here. Assuming the change in the phase ϕ s to be slower than the swimmer activity, we can use method of averaging. After substituting the aforementioned expressions for x i and y i in Eq. ( 25) and taking average over one cycle of the active force, we obtain following Adler equation

This demonstrates that for different α and β (which correspond to the distribution of the active force on the amoeboid swimmer surface) the phase-locking is observed under shear flow. The corresponding window of the the shear rates, however, is dependent on the nature of the active force distribution (α and β). It can be seen that for α = β = 0 (no propulsion) we get A c = A s = 0, implying lack of phase-locking.

FIG. 1 .

 1 FIG.1. Schematics of bead-spring and amoeboid microswimmers in shear flow.

FIG. 2 .

 2 FIG. 2. (A) Intrinsic viscosities of passive vesicles (f = 0) and amoeboid microswimmers (f > 0) for Γ = 0.75. γc marks the shear rate at which ∆η = 0. Inset shows [η] for different Γ with f = 40. (B) For small f , the normalization of the shear rate difference γ/ γc -1 and the active contribution to viscosity ∆η/ [η] p by f results in a reasonable collapse of the curves, especially in the central region; λ = 20.

FIG. 3 .

 3 FIG. 3. (A)The intrinsic viscosities (|χ|/fa = 0, 0.04 without (solid) and with (dashed) noise, respectively) and (B) normalized active viscosities, which follow Eq. (15) in the phase-locked regime, for the bead-spring triangles.

FIG. 4 .FIG. 5 .

 45 FIG.4. Dependence of (A) ωs, the angular frequency of the shear flow driven dynamics and (B) the average stresslet Σ for the bead-spring swimmer, as a function of shear rate, γ

FIG. 6 .

 6 FIG.6. Phase-locked regimes of the amoeboid swimmers for different f . ∆η/ [η] p for f = 5.0 is shown as a reference.

FIG. 7 .

 7 FIG. 7. Probability density of the presence of a bead (for triangle, upper panels) and swimmer membrane (for amoeboid, lower panels) under shear flow in the reference frame co-translating with the swimmer center.

2µ 54 (v b f a k 2 l 2 2µ108 ( 2 +

 5422 cos(α -β) + cos α + cos β + 3) + 4 √ 3 (9cζ + c/ζ) (sin β -sin α + sin(α -β)) + 729ζ 2 (9ζ + 1/ζ) (9ζ + 4/ζ) (9 + c 2 )(23) to be the effective suspension viscosity in the absence of coupling between activity and shear flow driven dynamics. It can be seen that for purely translation microswimmer (β = 0), we getη a ≈ η p -φ cos α) + 729ζ 2 (9ζ + 1/ζ) (9ζ + 4/ζ) (9 + c 2 )(24)

  dδ dt ≈ -( γ -ω a ) + f a k γ 9ζ + 4/ζ (A c cos δ + A s sin δ)(27)withA c = -2(2 cos α -3ζ sin α) + 3 √ 3ζ(cos β -1) + 2(cos β + 1) + sin β(3ζ + 2 √ 3)(28)A s = -2(2 sin α + 3ζ cos α) + 2 √ 3ζ(cos β -1) + 3ζ(cos β + 1) + sin β(2 -3 √ 3ζ).

ACKNOWLEDGEMENTS

We thank CNES (Centre National d'Etudes Spatiales) and the French-German University Programme "Living Fluids" (Grant CFDA-Q1-14) for the financial support. The simulations were performed on the Cactus cluster of the CIMENT infrastructure supported by the Rhône-Alpes region (GRANT CPER07 13 CIRA).

APPENDIX A. Numerical simulation of amoeboid microswimmer using Boundary Integral Method

We consider a single swimmer in a linear shear flow

Considering very small Reynolds number, the flow can be regarded as a Stokes flow and we can use the boundary integral method (BIM) [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF] to convert the Stokes equations into an integral form over the swimmer shape. The velocity at any point r in the domain can then be calculated by the boundary integral equation

where G and T are the single and double layer Green's functions, respectively, and the integrals are performed over the swimmer surface. The term Λ is defined as

Once the velocity on the membrane is known, the evolution of vesicle shape is obtained from a simple fixed time step Euler scheme

where u k (t) is obtained by taking Fourier transform of the velocity calculated from Eq. ( 18) at time t. Ideally, fluid incompressibility and membrane impermeability should automatically ensure the swimmer volume to be a constant but a drift due to numerical scheme can not be ruled out. To correct this drift, we inflate or deflate the swimmer through homogeneous normal deformation. We performed several confirmatory simulations with refined meshes and time steps, sampling points and Fourier harmonics to ensure the numerical accuracy of the scheme. This verification led us to compromise between efficiency and accuracy and characterize each swimmer surface by 63 Fourier harmonics. Further, we used 2048 sampling points to resolve the short-range hydrodynamics interactions. We calculated the swimmer surface velocity at 128 sampling points.