J Maysonnave

K Maussang
email: kenneth.maussang@lpa.ens.fr

N Jukam

J R Freeman

J Madéo

P Cavalié

R Rungsawang

S P Khanna

E H Linfield

A G Davies

H E Beere

D A Ritchie

S S Dhillon

J Tignon

Mode-locking of a terahertz laser by direct phase synchronization

Keywords: Semiconductor lasers, quantum cascade; (300.6495) Spectroscopy : Spectroscopy, teraherz

A novel scheme to achieve mode-locking of a multimode laser is demonstrated. Traditional methods to produce ultrashort laser pulses are based on modulating the cavity gain or losses at the cavity roundtrip frequency, favoring the pulsed emission. Here, we rather directly act on the phases of the modes, resulting in constructive interference for the appropriated phase relationship. This was performed on a terahertz quantum cascade laser by multimode injection seeding with an external terahertz pulse, resulting in phase mode-locked terahertz laser pulses of 9ps duration, characterized unambiguously in the time domain.

Introduction

Ultrashort pulses of light have attracted great interest for decades [START_REF] Hargrove | Locking of He-Ne laser modes induces by synchronous intracavity modulation[END_REF][START_REF] Ippen | Passive mode locking of the cw dye laser[END_REF][START_REF] Demaria | Self mode-locking of lasers with saturable absorbers[END_REF][START_REF] Keller | Recent Developments in compact ultrafast lasers[END_REF], resulting in major fundamental and technological advances, from metrology [START_REF] Udem | Optical frequency metrology[END_REF] to industrial applications [START_REF] Keller | Recent Developments in compact ultrafast lasers[END_REF]. They result from the temporal interference between modes that are equally spaced in the frequency domain by the free spectral range (FSR) of a cavity (Fig. 1a). For N modes with a random phase relationship, the combined intensity is random as well, with mean amplitude scaling as N (Fig. 1b, green curve). If, however, the phase difference between consecutive modes is fixed and constant in time, the laser is described as modelocked. The interference is then constructive, leading to the formation of a pulse train with peak intensity scaling as N 2 and period T given by the round trip time (Fig. 1b, blue curve). The experimental realization of mode-locking is based on two key elements [START_REF] Brabec | Intense few-cycle laser fields: Frontiers of nonlinear optics[END_REF]: first, a broadband gain medium that provides a large number of lasing modes leading to shorter pulses (~T/N); second, a mechanism that fixes the phase relationship between lasing modes. This is usually performed through internal modulation of the cavity, either active (e.g. using an acousto-optic modulator) or passive (e.g. the optical Kerr effect). Additionally, the fixed phase relationship may also be disrupted by the group velocity dispersion (GVD) of the medium. This can be compensated to achieve mode-locking by, for example, using an anomalous dispersive element.

Since the original demonstration of mode-locking [START_REF] Hargrove | Locking of He-Ne laser modes induces by synchronous intracavity modulation[END_REF], there has been a constant effort to realize similar schemes in various regions of the electromagnetic spectrum. The terahertz (THz) frequency domain, corresponding to energies in the milli-electronvolt range, has enormous potential for the exploration of intersubband physics in semiconductors, superconductors, imaging, gas analysis or biological applications [START_REF] Mittleman | Sensing with Terahertz Radiation[END_REF][START_REF] Tonouchi | Cutting-edge terahertz technology[END_REF]. In 2002, quantum cascade lasers (QCL) emerge as compact, powerful and efficient sources of THz radiation [START_REF] Köhler | Terahertz semiconductor-heterostructure laser[END_REF]. Over the past decade, QCL performance has steadily increased to give a far greater range of operating frequencies and higher operating temperatures. More recently, advanced functionalities have been investigated, including pulse amplification [START_REF] Jukam | Terahertz amplifier based on gain switching in a quantum cascade laser[END_REF], phase-control [START_REF] Oustinov | Phase seeding of a quantum cascade laser[END_REF][START_REF] Barbieri | Phase-locking of a 2.7 THz quantum cascade laser to a mode-locked erbium-doped fibre laser[END_REF], and mode-locking [START_REF] Barbieri | Coherent sampling of active mode-locked terahertz quantum cascade lasers and frequency synthesis[END_REF][START_REF] Wang | Mode-locked pulses from mid-infrared quantum cascade lasers[END_REF] offering the possibility of using QCLs as a powerful and versatile source for time-domain spectroscopy (TDS). Only one realization of mode-locking has been reported in the terahertz spectral range, with an active scheme based on frequency feedback [START_REF] Barbieri | Coherent sampling of active mode-locked terahertz quantum cascade lasers and frequency synthesis[END_REF] but, without direct control over the individual phases.

In the present work, thanks to the unprecedented control over the phase in the THz domain, we demonstrate a novel approach for mode-locking, where the phase of the laser longitudinal modes is directly imprinted using multimode injection seeding [START_REF] Oustinov | Phase seeding of a quantum cascade laser[END_REF]. We apply this technique to a 2.75THz quantum cascade laser (QCL), leading to the formation of 9ps THz laser pulses. The unique ability of Time Domain Spectroscopy [START_REF] Schmuttenmaer | Exploring Dynamics in the Far-Infrared with Terahertz Spectroscopy[END_REF] (TDS) and coherent electro-optic detection [START_REF] Kröll | Phase-resolved measurements of stimulated emission in a laser[END_REF][START_REF] Wu | Free-space electro optic sampling of terahertz beams[END_REF] to give access to both the amplitude and phase allows to fully characterize the locked phase relationship between the modes, as well as the initial transient laser dynamics. Corresponding time domain intensity profile. A given amplitude spectrum results in different behaviors in the time domain depending on the phase relationship between modes. The case of mode-locking, where all modes have the same phase, is illustrated in blue, whereas the case of a random phase distribution is shown in green.

Phase synchronization by multimode injection seeding

The principle of our experiment is schematically depicted in Fig. 2 and satisfies the two key requirements cited above. First, we use a QCL with a 350GHz gain bandwidth, able to sustain about 11 lasing modes. Second, the phase of the QCL lasing modes is initially fixed by an external seed with a known phase rather than set randomly by the inherent spontaneous emission [START_REF] Oustinov | Phase seeding of a quantum cascade laser[END_REF]. For this, the QCL is switched-on periodically and synchronously injected with broadband phase-locked THz seed pulses that have a constant phase over the QCL gain bandwidth [START_REF] Schmuttenmaer | Exploring Dynamics in the Far-Infrared with Terahertz Spectroscopy[END_REF]. This multimode injection seeding induces the formation of multiple longitudinal modes lasing with the fixed phase relationship of the seed. As a result, after a brief linear amplification regime, the gain is clamped to the losses, leading to the formation of a train of THz laser pulses with constant amplitude (Fig. 2a). Moreover, the QCL is periodically switched on and off, on a nanosecond timescale so that, as it will be demonstrated below, the dispersion remains negligible. This periodic switching is performed synchronously with the seeding, leading to a series of pulse trains, synchronized and phase-locked to each other. Fig. 2. a) Time domain: the QCL is gated at low frequency compared to the roundtrip time. When the laser is switched on, a seeding pulse is injected in the cavity with a field Eseed(t). After a brief transient regime, the laser gain g(t) is saturated and matches the losses, leading to an output field Eout(t) consisting in a constant amplitude laser pulse train. b) Frequency domain: the terahertz seed pulses consist in a frequency comb with spacing frep and no frequency offset (green full lines). The free running QCL modes have linewidths smaller than frep (blue full lines) and a priori do not match the excitation comb frequencies. When the laser is gated for Tgate< 1/frep, lines are broaden up to 1/Tgate with overlap over several frequencies of the excitation comb (blue dashed lines). If the laser is gated periodically at frep frequency, an inner comb structure appears, with frep spacing and no frequency offset (red full lines). Each component of the periodically gated QCL spectrum is then excited by a component of the excitation terahertz frequency comb, imprinting its phase.

The periodic switching also plays an important role for the resonant injection of the QCL laser modes, as can be understood in the frequency domain (Fig. 2b). The seed field consists in a series of identical THz pulses with a repetition rate f rep (76MHz). As a result, the seed spectrum is a frequency comb, with spacing f rep (Fig. 2c , green lines) and no frequency offset since the carrier phase shift is null [START_REF] Udem | Optical frequency metrology[END_REF]. The free-running QCL spectrum displays several modes, spaced by the FSR, with a linewidth in the MHz range [START_REF] Barkan | Linewidth and tuning characteristics of terahertz quantum cascade lasers[END_REF] (blue solid lines). Consequently, the seed pulse frequency components are a priori not resonant with the freerunning QCL modes. However, when the QCL is turned-on for a time T gate <1/f rep , each QCL mode broadens to a linewidth T gate (blue dashed lines), larger than the excitation comb frequency spacing. Last, since the QCL is periodically turned-on with a gate synchronized to seeding pulses, the spectrum of the broadened QCL is structured into a comb with spacing f rep and no frequency offset (red lines). Therefore each component of the periodically gated QCL spectrum is resonantly excited by a component of the seeding terahertz frequency comb, imprinting its own phase and permitting multi-mode injection seeding.

Experimental setup

A schematic diagram of the experiment is shown in Fig. 3. The QCL is an AlGaAs/GaAs bound-to-continuum design with a plasmon waveguide that lases at 2.75THz. The QCL is cleaved to a length of 1.5mm and soldered with indium onto a copper heat sink which is maintained at a temperature of 10K. An upper limit of the GVD was estimated to be 4210 -2 ps 2 /mm at 3THz. Thus, a gaussian pulse of 10ps FWHM broadens by one cycle (333fs) only after 3ns of propagation. Consequently, the QCL is switched-on for less than 1ns to avoid any dispersion effects. Fig. 3. a) The seed THz pulses are produced by illuminating an interdigited photo-conductive antenna with near infrared femtosecond laser pulses. After a THz seed pulse enters the QCL, a square RF pulse synchronized to the femtosecond laser pulse biases the QCL above lasing threshold. The THz beam emitted and a split-off beam from the femtosecond beam are then focused onto a ZnTe crystal for the electro-optic sampling detection. b) Detection of the emitted field, with a quarter wave electro-optic sampling geometry. The signal is proportional to the amplitude of THz field emitted by the QCL.

A square shape RF is generated by illuminating a photodiode with a Ti:Sapphire femtosecond laser, with pulse duration of typically 100fs and a repetition rate of 76MHz. This periodically pulses biases the QCL above lasing threshold [START_REF] Oustinov | Phase seeding of a quantum cascade laser[END_REF], when added to a quasi-DC bias (25% duty cycle at kHz frequencies) with a bias tee. The RF pulse rise time is typically 400ps with duration of 600ps. A 50Ω RF switch modulates this bias for lock-in detection. The THz seed pulses are generated with a photo-conductive antenna [START_REF] Madéo | Frequency tunable terahertz interdigitated photoconductive antennas[END_REF][START_REF] Dreyhaupt | High-intensity terahertz radiation from a microstructured large-area photoconductor[END_REF], excited by the same laser, and are coupled into the facet of the QCL with parabolic mirrors (not shown). An electronic delay on the RF pulse permits synchronization with the seeding pulse. The antenna is biased periodically at frequency 2f with a constant polarity. The QCL is modulated at frequency f=10kHz, whilst the lock-in reference is f, so that in this configuration only the field emitted by the QCL is detected. The THz beam emitted and a split-off beam from the femtosecond beam are then focused onto a 200 µm ZnTe crystal for the electro-optic sampling detection [START_REF] Wu | Free-space electro optic sampling of terahertz beams[END_REF] (Fig. 3b): the THz electric field induces a small birefringence and a polarization rotation of the femtosecond beam, measured in a balanced photodiode scheme. It results in a direct measurement of the field amplitude and phase, including any transient dynamics.

Results

Fig. 4. a) Amplitude and phase resolved measurements of the emitted electric field, for a THz seed pulse generated with 4V applied to the photo-conductive antenna, and a DC bias of 120A/cm 2 and 7.5W RF pulse applied to the QCL. The seeding pulse is injected in the QCL at t=-24ps. Inset: expanded view of the TDS measurement between t=526ps and t=528ps. b) Instantaneous emitted intensity, calculated from the TDS field measurement. A train of pulses, separated by the cavity round trip (35ps), is clearly visible, with a FWHM of 9ps in intensity. c) Pulse width (blue square, right axis) and spectrum width (red circle, left axis) as a function of the number of round trip in the QCL. d) Fourier transform of the measured electrical field in the steady state regime. Several modes are visible, separated by the cavity free spectral range (29GHz). A spectrum of the QCL, obtained by Fourier transform infrared spectroscopy (FTIR), has been superimposed. Both are renormalized to their maximum value and are resolution limited.

Figure 4a shows the amplitude of electrical field E(t) emitted by the QCL, and demonstrates unambiguously the formation of ultrashort THz laser pulses. At early times (t<300ps), the frequency components of the seed pulse corresponding to the QCL gain bandwidth are amplified, whilst the rest undergo losses and rapidly vanish, resulting in a spectral narrowing (blue squares in Fig. 4c). In the time domain, this results in pulse width broadening (red circles in Fig. 4c). After this transient regime of about 300ps, a quasi-constant amplitude train of pulses is clearly visible, with a FWHM of about 9ps in intensity, as illustrated more clearly by Fig. 4b which shows the square of the electric field E 2 (t). Figure 4d shows the spectrum obtained by Fourier transform of the signal measured in the steady state region (red curve). It is obtained by weighting the signal with a hamming window of 500ps centered at t 0 =550ps, and the FFT routine of MATLAB is used to compute the spectrum. It consists of an ensemble of discrete modes, equally spaced by the QCL FSR (29GHz). The resolution of the measurement is limited by the length of the mechanical delay line (800ps) that does not allow the resolution of the underlying QCL frequency comb discussed above. The spectrum of the free-running QCL (i.e. without seeding and RF modulation), which cannot be determined using electro-optic sampling, was measured with a commercial FTIR and is also shown in Fig. 4d (resolution limited). Whilst the center frequencies have similar relative amplitudes, the spectrum of the injected QCL presents additional lasing side-modes, interpreted as a signature of mode-locking [START_REF] Wang | Mode-locked pulses from mid-infrared quantum cascade lasers[END_REF]. Fig. 5. a) The phase of each laser mode (in blue) is obtained from the argument of the complex FFT of the field measured in the steady state region (t>300ps). Similarly, a measurement of the seeding pulse field provide after FFT the spectral phase of the seeding pulse (green line). In the inset, an expanded view for clarity, around the laser modes' frequencies, showing unambiguously that the phases are identical. The spectral phase of the seeding pulse has negligible variation over the frequency bandwidth of the laser (2π×6210 -3 rad standard deviation). b) Reconstructed field from the measured amplitudes and phases.

The spectral phase of the field can be obtained from the argument of the complex spectrum at the corresponding frequency, but for a window of 400ps only, centered at t 0 =550ps. Error bars are estimated by measuring the phase with the same procedure but for different values of t 0 from 500ps to 600ps. Then, the phase of each mode can be extracted (Fig. 5a, blue dots) and compared to the spectral phase of the seeding pulse (Fig. 5a, green line). They are measured to be almost identical within the error bars. This clearly demonstrates that (i) there is a fixed phase relationship between all the modes, imprinted by the seeding, and (ii) the phase difference between the modes is constant, providing optimal contrast for the interference and hence Fourier-transform limited pulses [START_REF] Brabec | Intense few-cycle laser fields: Frontiers of nonlinear optics[END_REF]. As a last demonstration of the completeness of the measurement, we show in Fig. 5b the field calculated from the measured amplitudes and phases (t>300ps), showing excellent agreement with the measured field (Fig. 4a).

Conclusion and outlook

In this current work, 9ps THz pulses are generated, limited by the number of modes i.e. the spectral gain of the laser. This can be easily increased using broadband QCLs [START_REF] Tur1inková | Ultra-broadband heterogeneous quantum cascade laser emitting from 2.2 to 3.2 THz[END_REF] that have shown bandwidths greater than 1THz. Further, new THz phase-modulation techniques [START_REF] Chen | A metamaterial solidstate terahertz phase modulator[END_REF] could be used to finely adjust the individual phases of the mode phases with on-demand values. This would provide an integrated pulse shaper, whereas in optics, pulse shaping is a subsequent [START_REF] Warren | Coherent control of quantum dynamics -the dream is alive[END_REF]. In conclusion, we have shown that mode-locking can be achieved by directly imprinting the desired phase relationship onto the laser modes of a THz QCL by multimode phase-locked injection seeding, directly characterized with phase resolved electro-optic detection. This could lead to experimental investigations of pulse formation dynamics and stability in QCLs [START_REF] Wang | Coherent instabilities in a semiconductor laser with fast gain recovery[END_REF][START_REF] Gordon | Multimode regimes in quantum cascade lasers: From coherent instabilities to spatial hole burning[END_REF][START_REF] Choi | Gain recovery dynamics and photon-driven transport in Cascade lasers[END_REF][START_REF] Gkortsas | Dynamics of actively mode-locked quantum cascade lasers[END_REF], and ultimately the generation of ultrashort and intense THz pulses. Such devices would also be powerful tools for use in THz-TDS, molecular spectroscopy or coherent control with phase-shaped pulses [START_REF] Warren | Coherent control of quantum dynamics -the dream is alive[END_REF][START_REF] Silberberg | Quantum Coherent Control for Nonlinear Spectroscopy and Microscopy[END_REF].

Fig. 1 .

 1 Fig. 1. a) Frequency domain representation for N=12 modes with equal amplitude. b)Corresponding time domain intensity profile. A given amplitude spectrum results in different behaviors in the time domain depending on the phase relationship between modes. The case of mode-locking, where all modes have the same phase, is illustrated in blue, whereas the case of a random phase distribution is shown in green.

Acknowledgments

We thank Louis-Anne de Vaulchier and François-Régis Jasnot, who performed the FTIR measurement. We acknowledge funding from ANR (FR), EPSRC (UK) and European Research Council grants "TOSCA" and "NOTES".


Julien Madéo

Pierrick Cavalié

Joshua R Freeman

Nathan Jukam

Jean Maysonnave

Kenneth Maussang

Harvey E Beere

David A Ritchie

Carlo Sirtori

Jérôme Tignon

Sukhdeep S Dhillon

All-optical wavelength shifting in a semiconductor laser using resonant nonlinearities

 Here, we demonstrate all-optical wavelength conversion of a near-infrared beam using a resonant non-linear process within a terahertz (THz) quantum-cascade-laser (QCL). 4 The process is based on injecting a low power CW near-infrared beam in resonance with the interband transitions of the QCL. This results in an enhanced nonlinearity allowing the efficient generation of the difference and sum frequency, shifting the frequency of the near-infrared beam by the QCL frequency. Efficiencies of 0.13% are shown which are equivalent to those obtained using Free Electron Lasers. As well as important implications as an ultrafast wavelength shifter, this work also opens up the possibility of efficiently up-converting THz radiation to the near-infrared and the study of high THzoptical field interactions with quantum structures using QCLs.

Wavelength division multiplexing (WDM) is currently used extensively in optical fibre networks to enhance the carrier capacity of optical fibres where each wavelength in a multiwavelength bunch is assigned to a particular communication channel. [START_REF] Decusatis | Handbook of Fiber Optic data Communication[END_REF] In these types of networks, wavelength manipulation is essential for data routing, and optoelectronic shifters are employed to perform wavelength shifting; the optical signal is converted to an electrical signal and then back into an optical signal at a different wavelength. This creates an undesirable speed bottleneck. In order to overcome this problem, as well as to increase the bandwidth, all-optical networks have been proposed that would succeed their electrical counterparts. For these types of networks new types of optical devices are required that directly manipulate communication channels to shift their wavelength across the optical fiber bandwidth.

Non-linear processes in semiconductor devices have the potential to fill this technological gap where two wavelengths in a nonlinear material can be mixed to generate the sum or difference frequency, effectively shifting the original wavelengths. [START_REF] Campi | Wavelength conversion technologies[END_REF] Normally, however, these techniques are based on the small bulk non-linearities of the material, and demand phase matching for a long interaction length and/or high pump powers. These considerations can be overcome by using resonant non-linearities of quantum wells [START_REF] Sirtori | Far-infrared generation of doubly resonant difference frequency mixing in a coupled quantum well two dimensional electron gas system[END_REF][START_REF] Belkin | Terahertz quantum-cascade-laser source based on intracavity differencefrequency generation[END_REF] which are orders of magnitudes greater than those of the bulk, permitting much shorter interaction lengths. Indeed considerable investigations have been undertaken of the wavelength shift of near infrared (NIR) beam in the presence of an intense THz beam in a quantum well system. 8- 11 These studies were based on enhanced non-linear susceptibilities where the near-infrared beam was resonant with excitonic interband transitions and the THz beam was resonant with excitonic intersubband levels. The resonances could be modified to engineer the wavelength shift by changing the quantum well geometry [START_REF] Phillips | Generation of first-order terahertz optical sidebands in asymmetric coupled quantum wells[END_REF] or by the application of an external electric field [START_REF] Ciulin | Terahertz optical mixing in biased GaAs single quantum wells[END_REF][START_REF] Su | Voltagecontrolled wavelength conversion by terahertz electro-optic modulation in double quantum wells[END_REF] and large efficiencies (0.1% -0.2%) could be obtained. [START_REF] Wagner | Resonant enhancement of second order sideband generation for intraexcitonic transitions in GaAs/AlGaAs multiple quantum wells[END_REF][START_REF] Carter | Terahertz electro-optic wavelength conversion in GaAs quantum wells: Improved efficiency and room-temperature operation[END_REF] However, an important point of all this work is that the THz beam is provided by an entire facility -the free electron laser (FEL) -and thus prohibiting its relevance to real world applications. THz QCLs [START_REF] Köhler | Terahertz semiconductor-heterostructure laser[END_REF][START_REF] Williams | Terahertz quantum cascade lasers[END_REF] are recently realised semiconductor sources that operate in the THz range. QCLs are based on intersubband transitions where laser action takes place between confined conduction band subbands in a series of coupled quantum wells. The intracavity fields (upto a few kV/cm) of these devices can approach those that are used in the FEL studies mentioned above. [START_REF] Carter | Terahertz electro-optic wavelength conversion in GaAs quantum wells: Improved efficiency and room-temperature operation[END_REF] QCLs therefore have the potential as an integrated wavelength converter to shift an external NIR beam by the QCL frequency by (i) providing large intracavity powers and (ii) an enhanced non-linearity from the interband resonance of the NIR beam with the confined states. Previous work using QCLs and a NIR beam took advantage of the bulk second order non-linearity of GaAs to perform frequency mixing [START_REF] Dhillon | Terahertz transfer onto a telecom optical carrier[END_REF] or the use of a double resonant process [START_REF] Zervos | Coherent near-infrared wavelength conversion in semiconductor quantum cascade lasers[END_REF] but the efficiencies of the processes were limited (10 -4 -10 -3 %). Here, we demonstrate an improvement of 2 orders of magnitude with conversion efficiencies of up to 0.13 %, comparable to those obtained in FEL investigations.

Figure 1a shows the schematic of the process investigated here via the resonant interband excitation of the QCL with the NIR beam. The THz QCL laser transition E QCL occurs within the conduction band between the highlighted green states (green wave arrow).

A NIR beam E NIR (red arrow) is coupled into the QCL cavity and resonantly tuned with an interband transition implying hole and electron states. As a result the difference frequency is generated E NIR -E QCL (dark red arrow), via a virtual state and the THz photon (green arrow), which is below the bandgap and therefore avoids absorption. [START_REF] Phillips | Generation of first-order terahertz optical sidebands in asymmetric coupled quantum wells[END_REF] (The reverse situation also occurs with an excitation at the virtual transition to generate the sum frequency, E NIR +E QCL , at the bandgap).

Figure 1b shows the geometry of the experiment where the input interband excitation and the THz QCL emission are collinear i.e. in the same plane parallel to the surface of the QCL. This is in strong contrast to previous experiments that have investigated resonant nonlinear mixing where the THz beam and near-infrared excitation were orthogonal. [START_REF] Phillips | Generation of first-order terahertz optical sidebands in asymmetric coupled quantum wells[END_REF][START_REF] Ciulin | Terahertz optical mixing in biased GaAs single quantum wells[END_REF][START_REF] Carter | Terahertz electro-optic wavelength conversion in GaAs quantum wells: Improved efficiency and room-temperature operation[END_REF][START_REF] Zervos | Coherent near-infrared wavelength conversion in semiconductor quantum cascade lasers[END_REF] This type of guided geometry for both the THz beam and input interband excitation allows the use of a much longer interaction length. The NIR pump is coupled into one QCL facet and the difference frequency exits the opposite facet as well as the remaining input (the latter depending on the optical losses).

The guided modes are shown in Figure 1c for the THz QCL emission and the injected NIR beam (just below bandgap) using a dual wavelength waveguide. The confinement of the transverse magnetic (TM) polarised THz beam is based on a standard surface plasmon mode.

The NIR excitation (transverse electric (TE) or TM) is confined by the top metal layer and a 300nm Al 0.5 Ga 0.5 As grown between the lower doped layer and the substrate. This layer has a lower refractive index than the surrounding material and therefore dielectrically confines the injected NIR beam.

QCLs operating at 2.8THz using GaAs/AlGaAs quantum wells and based on a boundto-continuum design [START_REF] Barbieri | 9 THz quantum cascade lasers operating up to 70 K in continuous wave[END_REF] were employed and operated in continuous wave at 10K. (See method section for more details on the sample and experimental set-up). Figure 2a show the spectra of the beam without transmission through the QCL with the pump beam centred at E NIR =1.5267eV (λ=812nm) i.e. just above the effective bandgap of the QCL and corresponding to electronic transitions between the first confined hole and electron states. The polarisation of the NIR beam was chosen to be parallel to that of the QCL (i.e. TM polarized) implying interband transitions from only the light hole states. [START_REF] Rosencher | Optoelectronics[END_REF] Figure 2b shows the spectrum after transmission through the QCL driven below laser threshold (black curve) where E NIR is just visible, as a result of a parasitic part of the beam that does not pass through the ridge. The situation changes drastically when the QCL is above threshold (red curve). The difference frequency is clearly observed at E d =E NIR -E QCL = 1.5152eV (λ=818nm), i.e. separated from the pump E NIR by exactly the photon energy of the THz QCL (f=2.78THz) and is below the bandgap of the material. A point to note is the high apparent conversion efficiency, with the difference frequency 14 times more intense than the pump wavelength (taken as the ratio of the integrated signals of the pump and sideband). This is due to the sharp interband absorption. To estimate the actual conversion efficiency, defined as the ratio of the power of the sideband divided by that of the input NIR pump P d /P NIR , the coupling efficiency of the pump needs to be determined. This was done by characterising the transmission of the pump at an energy below the effective gap of the QCL where the interband losses are zero. This allows the determination of the NIR pump intensity coupled into the input facet of the QCL.

With this calibration taken into account an efficiency of 0.13% is determined.

The inset of figure 3a shows the interest of this technique for THz detection. A high resolution spectrum of the sideband (with the pump at 1.528eV) showing many modes is shown which is an exact replica of the QCL Fabry Perot emission and which was taken in less than a second. This shows that this wavelength conversion allows one to measure and upconvert the THz emission of the QCL to the near-infrared. Further, these measurements used standard CCD camera technology for detection. The spectrum exhibits a 1.2 GHz resolution which is comparable to that of a considerably slower high resolution FTIR spectrometer.

The resonant nature of interaction can be seen in figure 3a where spectra for several pump wavelengths and their corresponding difference frequency for a TM pump polarisation are shown. For clarity, the curves have been normalised by setting the pump wavelength intensities to one. As the pump energy is increased from 1.522eV to 1.534eV the difference frequency increases in intensity showing a double resonance before decreasing at higher pump energies. These results are plotted in figure 3b and show the efficiency with the absorption of the pump taken into account (see above) as a function of the pump energy (square points). Resonances at 1.531eV and 1.527eV are observed with the latter showing a conversion efficiency of 0.13% which was the highest obtained. This is more than two orders of magnitude greater than previously demonstrated frequency mixing using QCLs. The full width at half maximum of each resonance is ~1 -2 meV. These resonances arise when E NIR is resonant with interband transitions that have a large overlap between the electron and hole wavefunctions. Also shown in figure 3b is the efficiency with a TE polarised NIR beam corresponding to interband excitation predominately from heavy holes states. [START_REF] Rosencher | Optoelectronics[END_REF] One resonance is clearly observed at 1.525 eV and another broader one around 1.528 eV. The slightly lower energy of the resonance for the TE polarization compared to that of the TM is due to the higher confinement of the heavy hole states. The slighty reduced efficiency in the TE polarisation is due to the smaller overlap of the electron and heavy hole states.

To identify the states involved in the nonlinear process, the photoluminescence (PL) from the QCL was investigated (Figure 4b). As the PL emission is given by the transitions having the lowest energy, i.e. involving heavy holes states, it is compared to the efficiency curve for a TE polarization of the NIR pump (figure 4a, red squares). The TE efficiency is also plotted as a function of the difference energy E NIR -E QCL (green circles). Firstly, regarding the resonances in the efficiency at the pump beam energies (red squares), the main peak corresponds to the shoulder in the PL spectrum at 1.525 eV and that the higher energy and broader resonance corresponds to a transition around 1.527 eV. Secondly, no peaks are seen in the PL in the difference energy range (green circles) illustrating that there are no resonant transitions and confirming that the wavelength conversion is realized through a virtual state below the bandgap. (The small shoulder at 1.519eV is due to PL from the GaAs substrate).

Figure 4c shows a comparison between the PL spectrum with the overlap of the interband wavefunctions of the QCL, the PL being proportional to the square of the overlap between the electron and hole states. [START_REF] Freeman | Photo-luminescence study of heterogeneous terahertz quantum cascade lasers[END_REF][START_REF] Vitiello | Measurement of subband electronic temperatures and population inversion in THz quantum cascade lasers[END_REF] Figure 4d presents the bandstructure of the studied sample including conduction and valence bands. The states taken into account are highlighted and are the three lowest lying heavy hole levels (labelled H1 to H3) and five electronic states in the miniband (labelled E1 to E5) which have a significant overlap. Comparing figures 4a, 4b and 4c, the efficiency peak at 1.525eV is a result of resonances with HH1E1 and HH1E2 transitions i.e. the lowest lying states. The higher energy broader resonance is more difficult to identify to a sole transition and appears to be related multiple contributions between the three hole states and the lower electron states. This explains the broader nature of the second resonance. A similar analysis can be performed for the light hole states that shows only a significant overlap comes between the lowest lying hole state and the lowest electronic states that are separated by 3meV, in agreement with the results of figure 3b. Thus the nonlinear process is singly resonant with the lowest energy hole states and the electronic states of the QCL's miniband.

It is possible to estimate the second order susceptibility, χ [START_REF] Desurvire | Capacity demand and Technology Challenges for Lightwave Systems in the Next Two Decades[END_REF] , from the efficiency, η: where P d (n d ), P NIR (n NIR ), P QCL (n QCL ) are the intensities (refractive indices (~ 3.6)) of the generated beam, the input NIR pump, and the THz QCL respectively, L is the cavity length of the QCL (1.5mm), λ d is the generated wavelength (818nm) and S is the interaction area defined as the modal overlap between the three interacting waves (8000 µm 2 ) [17]. The intracavity THz power is 75mW, estimated from the detected output power and a facet reflectivity of 0.32. p are the losses of the NIR pump (estimated at 1000 cm -1 ) [START_REF] Rosencher | Optoelectronics[END_REF] and k is the phase mismatch. Here the losses of the difference frequency and the THz beam are taken as zero. Although phase matching in the geometry presented here is possible due to refractive index dispersion at resonance [START_REF] Blakemore | Semiconducting and other major properties of gallium arsenide[END_REF] (see supplementary material S2), the high losses of the pump beam is the limiting factor regarding the interaction length. Equation 1 can then be simplified to:

23- 25 2 2 2 2 2 2 0 2 2 ) 2 ( 2 4 2 4 2 sin 8       +       Δ       +       Δ = = - L kL L sh kL e S c
p d d QCL NIR QCL NIR d S c n n n P P P α λ ε χ π η 1 64 2 0 2 ) 2 ( 2 ≈ = [2]
Taking the maximum efficiency of 0.13 %, a second order non-linearity of ~1×10 4 pm/V is determined. This is rough agreement with previous studies which have shown interband nonlinearities in the range of 10 2 -10 4 pm/V for quantum wells [START_REF] Garmire | Nonlinear optics in semiconductors II[END_REF][START_REF] Khurgin | Second-order nonlinear effects in asymmetric quantum-well structures[END_REF] confirming that the nonlinear susceptibility is enhanced by the resonant excitation. [START_REF]This rather large range in simulated and experimental nonlinearities is a result of interband nonlinearities being a complex problem that involves the contribution of many real and virtual band-to-band transitions[END_REF] We have shown that the resonant interband properties of a QCL can be used for efficient frequency mixing through enhanced nonlinearities. The perspectives on this work are wide ranging. The wavelength shift can be engineered to any desired THz value [START_REF] Scalari | THz and sub-THz quantum cascade lasers[END_REF] and can be equally applied to MIR QCLs, [START_REF] Faist | Quantum Cascade Laser[END_REF] where the wavelength shift is much greater and can be used to shift between different telecommunication bands [START_REF] Shin | All fiber optical inter-band router for broadband wavelength division multiplexing[END_REF] . As well as room temperature and high power output, [START_REF] Lu | 4 W room temperature continuous wave operation of distributed feedback quantum cascade lasers[END_REF] MIR QCLs are based on InGaAs/AlInAs quantum wells where the interband transition is directly in the telecommunication range. Further increases in efficiencies could be realised through a) adapted active region designs or the insertion of passive quantum wells to enhance the non-linearity through optimisation of the overlaps between the confined states; and b) the combination of intersubband non-linearities with those of the interband transitions.

This work also shows the potential to efficiently up convert the QCL emission permitting the use of NIR technology for the detection of THz emission or providing a NIR-THz link for free space telecommunications. On a more fundamental side, this work also opens up the possibility of studying high THz-optical field interactions using compact and powerful QCLs, previously reserved to entire facilities such as the FEL.

In conclusion an efficient wavelength converter based on a compact QCL was demonstrated.

Frequency conversion with high efficiencies was realised through the enhancement of the non-linearity where the pump is resonant with interband transitions, combined with the high THz intra-cavity power density. These developments show the potential of QCLs as novel optical components for the all-optical telecommunication networks.

Methods

QCLs operating at 2.8THz using GaAs/AlGaAs quantum wells and based on a boundto-continuum design were employed (12µm active region thickness). Samples were processed into a single plasmon geometry with a ridge width of 250µm and a cavity length of 1.5mm.

The samples were operated in continuous wave at 10K. (See supplementary material S1 for the optical and electrical characteristics as well as the spectrum). As the bandgap of GaAs is in the near-infrared, the interband pump was sourced from a CW Ti:Sapphire laser that also allowed a large wavelength tunability and permits the correct interband resonance excitation to be found. 100 µW of NIR pump power was used with a coupling efficiency of ~20%, resulting in ~20µW coupled into the QCL cavity. Low powers for the input NIR beam were used so as not to affect the QCL performance. This was verified by confirming that the threshold current observed in the change in differential resistance of the VI did not increase with the coupled NIR beam. The transmitted NIR beam is collected at the opposite facet using a high numerical aperture objective and analysed using a spectrometer coupled to a thermoelectrically cooled CCD camera. Frequency (THz)

Figure Captions

Figure 1|

 1| Figure 1| Scheme principle and optical modes. a, Schematic of the resonant non-linear

Figure 2|

 2| Figure 2| Wavelength shifting using interband excitation. a, Spectrum of the NIR pump,

Figure

  Figure 3| Resonant behavior and polarization effect of NIR pump. a, Spectra of the QCL

Figure 4|

 4| Figure 4| Confined states involved in resonant nonlinear interaction. a, Conversion

E NIR -E QCL

					Frequency (THz)
	u.) Intensity (a.u	14 16 12 2 4 6 8 10	366 369 Frequency (THz) 367 367 368 370 ab 368 369 370 a 0.15 TM 368 369 370 371 Frequency (THz) 366 371 a. u.) Intensity (a TE %) 0.05 0.10 Efficiency (% 1.5164 1.5168 1.5172 Energy (eV) E g Pump
		0 0	1.515	1.520	1515 1.515 1.525 1.530	1520 1.520	1 525 1.525 1.520 0.00	1530 1.525 1.530	1.530	1.535
				E QCL ǻE= 11.5 meV ~ 2.78 THz Q E NIR -E QCL E NIR Pump Energy (eV) b Energy (eV) Intens sity (a.u.)
					1.515	1.520	1.525	1.530
					Energy (eV)
			Hole
			State
				Valence Band

Acknowledgements

This work was financially supported by the "Programme Francilien de Recherche en Nanosciences" (CNano-IDF, contract TeraConversion) and the "Agence Nationale de la Recherche" (ANR, contract "HI-TEQ" ANR-09-NANO-017). The author J.R.F would like to acknowledge funding from the Marie Curie Action fellowship (Grant no. 274602). The author J.M. would like to acknowledge funding from French Ministry of Defense (DGA). The "Laboratoire Pierre Aigrain" (LPA) is a "Unité Mixte de Recherche Associée à l'ENS, au Centre National de la Recherche Scientifique (CNRS) UMR8551 et aux Universités Paris 6 et 7". The device fabrication was performed at the nanocenter "La Centrale de Technologie Universitaire" (CTU-IEF-Minerve), which was partially funded by the "Conseil General de l'Essonne".


N Guizard

V S Fonov

B Aubert-Broche

D García-Lorenzo

P Coupé

S F Eskilden

D L Collins

Robust individual template pipeline processing for longitudinal MR images

Longitudinal measures of brain volume are powerful tools to assess the anatomical changes underlying on-going neurodegenerative processes. In different neurological disorders, such as in multiple sclerosis, Alzheimer's disease and Parkinson's disease, the neurodegenerative aspect may result in subtle anatomical brain changes before the appearance of clinical symptoms. Large longitudinal brain imaging datasets are now accessible to investigate such structural changes over time and to evaluate their use as biomarkers of prodromal disease. However, manual segmentation is long and tedious and although automatic methods exist, they are often performed in a cross-sectional manner where each visit is analysed independently. With such analysis methods, bias, error and longitudinal noise may be introduced. MR scanner noise and physiological effects can also introduce additional variability. Therefore, we developed a specific pipeline for longitudinal brain image analysis. To avoid any bias, an individual subject template is created and used as a reference within the pipeline. Then, the pair-wise deformation fields of each visit to the individual template are used to estimate the variation between individual time-points.

Introduction

Longitudinal measures of brain volumetry are powerful tools to assess the anatomical changes underlying on-going neurodegenerative processes. In different neurological disorders, such as multiple sclerosis (MS), Alzheimer's disease (AD) and Parkinson's disease (PD), brain atrophy has been shown to be a good surrogate marker of disease progression [START_REF] Chard | The longitudinal relation between brain lesion load and atrophy in multiple sclerosis: a 14 year follow up study[END_REF][START_REF] Burton | Cerebral atrophy in Parkinson's disease with and without dementia: a comparison with Alzheimer's disease, dementia with Lewy bodies and controls[END_REF][START_REF] Ridha | Tracking atrophy progression in familial Alzheimer's disease: a serial MRI study[END_REF]. Magnetic resonance imaging (MRI) can provide reproducible 3D structural images of the brain, which can be used to assess its integrity. Furthermore, the emergence of freely available longitudinal MRI databases, (e.g., ADNI [START_REF] Mueller | The Alzheimer's disease neuroimaging initiative[END_REF], OASIS [START_REF] Marcus | Open Access Series of Imaging Studies: Longitudinal MRI Data in Nondemented and Demented Older Adults[END_REF]) provide the necessary data to develop and test new methods and investigate the longitudinal structural changes of healthy and pathological brains.

MRI-based neuro-anatomical studies are often performed in a cross-sectional manner where each time-point is evaluated independently. Typically, brain morphometry comparison can be done by matching paired images (template-to-subject or subject-to-subject), where the deformation field is used to map atlas regions for ROI analysis, or for voxel-wise comparisons of anatomical changes such as deformation-based morphometry (DBM), or by segmenting anatomical structures from each individual visit independently. However, in the context of longitudinal datasets, the robust estimation of anatomical changes is still challenging [START_REF] Thompson | Bias in tensor based morphometry Stat-ROI measures may result in unrealistic power estimates[END_REF].

Different approaches have been proposed to overcome the complexity of anatomical longitudinal data image analysis. In the context of clinical evaluation over a few years where anatomical changes are small and continuous, the use of 3D individual template targets have been proposed to perform non-linear registration [START_REF] Reuter | Within-subject template estimation for unbiased longitudinal image analysis[END_REF]. Indeed, to compare anatomical differences, 3D population templates have proven their importance for different applications such as mapping (function, structure, vasculature, etc) [START_REF] Thompson | A framework for computational anatomy[END_REF], volume estimation [START_REF] Collins | ANIMAL: Validation and Applications of Non-Linear Registration-Based Segmentation[END_REF] and group comparisons [START_REF] Ashburner | Identifying global anatomical differences: Deformationbased morphometry[END_REF]. While different techniques exist to create unbiased population templates for cross-sectional studies [START_REF] Guimond | Automatic Computation of Average Brain Models[END_REF][START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF], few of these techniques have applied for the creation of an individual subject 3D template. Reuter et al. created a 3D template for longitudinal analysis by computing the median of the rigidly registered subject images [START_REF] Reuter | Within-subject template estimation for unbiased longitudinal image analysis[END_REF].

In this article, we propose a method to create an unbiased individual template from non-linear registration of the subject's visits. The method is evaluated using segmentations of the brain: lateral ventricles and hippocampi obtained by non-linear matching of the individual template to a population-specific template (i.e., an AD template), followed by a non-local patch-based segmentation technique. Inspired by the non-local means technique proposed by Buades et al. [START_REF] Buades | A non-local algorithm for image denoising[END_REF], the non-local patchbased segmentation uses the redundancy of similar regions from different subjects to identify and label corresponding local structure. This segmentation approach is inspired by the work of Coupe et al. [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF], Eskildsen et al. [START_REF] Eskildsen | BEaST: brain extraction based on nonlocal segmentation technique[END_REF] and Fonov et al [START_REF] Fonov | Automatic lateral ventricle segmentation in infant population with high risk of autism[END_REF].

Method

The pipeline developed in the context of longitudinal MRI segmentation consists of two parts (Fig. 1). In summary, the first steps of the pipeline are performed crosssectionally (CS) or independently for each of the subject visits. Then the preprocessed data goes through a subject-specific pipeline which considers the subject visits as an ensemble to generate the unbiased individual template for the subject i first using linear registration to build a subject-specific linear template (SL) and then using non-linear registration to build a subject-specific non-linear template (SNL). The method is detailed in the following sections.

2.1

Cross-sectional pre-processing MRI denoising: to improve the performance of the pipeline, the redundancy of the image is used to remove the noise introduce during the MRI acquisition. The standard deviation of the MRI Rician noise is estimated automatically following [START_REF] Coupé | Robust Rician noise estimation for MR images[END_REF] Intensity non-uniformity correction and normalization: a non-parametric estimation of the slow varying non-uniformity field corrects the intensity inhomogeneity produce by scanner radio-frequency coil variations [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF]. In addition, linear histogram matching is performed between each subject and a reference image (AD template) to normalize the image intensities between subjects/scans to a range between 0.0 and 100.0.

Stereotaxic space: to correct for variation in head position and orientation but to preserve the subject head size, an initial linear rigid registration (translation and rotation) is computed to bring each subject into the ICBM152 template stereotaxic space [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF].

Subject specific pipeline

Individual template creation: The subject-specific template is based on the work of Guimond et al. and Fonov et al. [START_REF] Guimond | Automatic Computation of Average Brain Models[END_REF][START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF][START_REF] Guimond | Three-dimensional multimodal brain warping using the demons algorithm and adaptive intensity corrections[END_REF]. The template is created in two steps, first using linear registration and second, using non-linear registration:

─ For each subject i, and each visit j, a hierarchical nine-parameter linear registration, based on an intensity cross-correlation similarity measure, is performed between the volume V i,j stx and the current linear template SL i [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF] and the first visit is used an initial target.

─ Removal of the average transformation for each subject transformation ( ) to remove bias ( ).

─ The new individual template is obtained by averaging each visit with their respective corrected transformations ( ).

─ To improve the inter-visit intensity inhomogeneity, at each iteration we compute the residual inhomogeneity field [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF] of the difference image between the resampled visit and the subject template.

After creation of SL i , the subject-specific linear template, the non-linear subjectspecific template SNL i is computed with a similar approach, but using non-linear registration, ANIMAL [START_REF] Collins | ANIMAL: Validation and Applications of Non-Linear Registration-Based Segmentation[END_REF] (Fig. 1). The template SL i is used as the target for the first iteration and the details of the hierarchical schedule is summarized in Table 1. Table 1. ANIMAL non-linear registration schedule. For each iterations, we define a step size as the distance between control nodes for the free-form deformation recovered. The blurring kernel is the size of the full-width-half-maximum of the Gaussian kernel used to blur the source and target data. The local correlation which define the local similarity is estimated in the neighboorhood of diameter equals to the neighboorhood size parameter.

Template segmentation and propagation:

The longitudinal pipeline applies the proposed methods of patch-based segmentation proposed initially by Coupé et al. [START_REF] Coupé | Patch-based segmentation using expert priors: application to hippocampus and ventricle segmentation[END_REF] for hippocampus. then proposed for ventricles by Fonov et al. [START_REF] Fonov | Automatic lateral ventricle segmentation in infant population with high risk of autism[END_REF] and the whole brain by Eskildsen et al. [START_REF] Eskildsen | BEaST: brain extraction based on nonlocal segmentation technique[END_REF]. In summery, based on the similarity of each patch of the SNL i and its surrounding patches, to all patches of the library in a certain search volume, a label is attributed to the considered voxel. Experts manually segmented the different regions considered, however different cohorts were chosen therefore each segmentation is performed independently.

Once the SNL i template is segmented, the non-linear deformations for each visit used to create the subject-specific template are then applied to transform the segmentations to each visit V i,j for the subject (Fig. 2).

Fig. 2. Lateral ventricles, hippocampi and brain patch-based segmentation. The subject nonlinear template is warpped into the ICBM template space to perform a non-local segmentation of the lateral ventricles, hippocampi and brain. The patch-based approach consists in finding the most similar voxels in a given neighbourhood of the T1W image library and then fuses the corresponding labels to create the segmentation of the different structures.

2.3

Data

The NIBAD challenge longitudinal data consists of T1W images (1.5x0.9375, 0.9375mm 3 , 1.5T MRI scanner) of 46 patients with AD and 23 age-match controls scanned at 0, 2, 6, 12, 26, 38 and 52 weeks.

Results

Without a gold standard provided for the NIBAD challenge, it is impossible to measure the accuracy and precision of the proposed pipeline. However, the patchbased technique for anatomical region segmentation has been shown to yield good results compare to the state of the art. Figure 3 is a random selection of 5 subject nonlinear templates with their corresponding brain, ventricles and hippocampi mask before propagation to the subject visits.

Discussion and conclusion

The proposed pipeline with individual template creation in the context of longitudinal analysis provides accurate segmentation of brain, ventricles and hippocampi. Despite the lack of temporal information, the pipeline treats the ensemble of subjects first independently before creating a robust individual 3D template for segmentation, and therefore avoiding non-symmetric and/or non-transitive errors.

The patch-based technique for segmentation presents the advantage of being accurate and also it can produce multiple segmentations at once, if the pre-labeled library used as references contains more than one labeled structure. While the segmentation library used here has been created from different cohorts of subjects (e.g., ADNI, ICBM, etc…), visual quality control of the NIBAD challenge data confirmed the accuracy of the segmentations.

In conclusion, our pipeline fulfills the requirements for the NIBAD challenge by providing brain, ventricles and hippocampi volumes and atrophies. All the longitudinal data completed the pipeline, including subjects with less conventional head positioning and poor quality acquisitions. If temporal information was provide, more sophisticated longitudinal regularization could improve the results by correcting for physiological, scanner physical variability, motion artifacts.

Fig. 1 .

 1 Fig. 1. Longitudinal pipeline. The different steps performed on each subject visits are represented in the left part of the diagram, where the processes in the left small square represents the cross-sectional (CS) part of the pipeline. The subject template creation (linear and non-linear) is represented in the right side of the figure.

Fig. 3 .

 3 Fig. 3. Individual non-linear template brain, ventricles and hippocampi segmentation for subjects 188, 192, 193, 230 and 249.


Nicolas Guizard

Vladimir S Fonov

Daniel García-Lorenzo

Bérengère Aubert-Broche

Simon Eskildsen

D Louis Collins

S F Eskilden

Spatio-temporal regularization for longitudinal registration to an unbiased 3D individual template

Keywords: Longitudinal registration, spatio-temporal consistency, unbiased template creation

Introduction

Longitudinal measures of brain volumetry are powerful tools to assess the anatomical changes underlying on-going neurodegenerative processes. In different neurological disorders, such as multiple sclerosis (MS), Alzheimer's disease (AD) and Parkinson's disease (PD), brain atrophy has been shown to be a good surrogate marker of disease progression [START_REF] Chard | The longitudinal relation between brain lesion load and atrophy in multiple sclerosis: a 14 year follow up study[END_REF][START_REF] Burton | Cerebral atrophy in Parkinson's disease with and without dementia: a comparison with Alzheimer's disease, dementia with Lewy bodies and controls[END_REF][START_REF] Ridha | Tracking atrophy progression in familial Alzheimer's disease: a serial MRI study[END_REF]. Magnetic resonance imaging (MRI) can provide reproducible 3D structural images of the brain, which can be used to assess its integrity. Furthermore, the emergence of freely available longitudinal MRI databases, (e.g., ADNI [START_REF] Mueller | The Alzheimer's disease neuroimaging initiative[END_REF], OASIS [START_REF] Marcus | Open Access Series of Imaging Studies: Longitudinal MRI Data in Nondemented and Demented Older Adults[END_REF]) provide the necessary data to develop and test new methods and investigate the longitudinal structural changes of healthy and pathological brains. MRI-based neuro-anatomical studies are often performed in a cross-sectional manner where each time-point is evaluated independently. Typically, brain morphometry comparison can be done by matching paired images (template-to-subject or subjectto-subject), where the deformation field is used for re-mapping atlas regions or for voxel-wise comparisons of anatomical changes such as deformation-based morphometry (DBM). However, in the context of longitudinal datasets, the robust estimation of anatomical changes is still challenging [START_REF] Thompson | Bias in tensor based morphometry Stat-ROI measures may result in unrealistic power estimates[END_REF]. Indeed, if we assume that longitudinal changes are smoothly varying, spatially local, monotonic temporal processes, consideration of individual time-points independently can generate unneeded noisy longitudinal measurements due to the intrinsic noise associated with each visit. Different studies have shown the impact of MRI acquisition on structural measurements [START_REF] Caramanos | Gradient distortions in MRI: Characterizing and correcting for their effects on SIENA-generated measures of brain volume change[END_REF] and cortical thickness [START_REF] Westlye | Increased sensitivity to effects of normal aging and Alzheimer's disease on cortical thickness by adjustment for local variability in gray/white contrast: A multi-sample MRI study[END_REF]. Therefore, methods, which integrate constraints from the temporal dimension (i.e., 4D methods), should produce more accurate, robust and stable measures of the longitudinal anatomical changes resulting in a more realistic estimation of temporal evolution.

Different approaches have been proposed to overcome the complexity of anatomical 4D longitudinal data image analysis. In the context of clinical evaluation over a few years where anatomical changes are small and continuous, the use of 3D individual template targets have been proposed to perform non-linear registration [START_REF] Reuter | Within-subject template estimation for unbiased longitudinal image analysis[END_REF].

More sophisticated and mathematically proven approaches have been proposed in the context of larger anatomical changes over time (i.e. over the span of childhood). For example, a 4D population model creation using Gaussian kernel regression has been suggested by Davis et al. [START_REF] Davis | Population Shape Regression from Random Design Data[END_REF] where each image is registered independently to a moving average, avoiding creating explicit parameterized mode of the longitudinal changes. Kernel regression has also been used in the framework of the Large Deformation Diffeomorphic Metric Mapping (LDDMM) [START_REF] Hart | DTI Longitudinal Atlas Construction as an Average of Growth Models[END_REF][START_REF] Prastawa | Building spatiotemporal anatomical models using joint 4-D segmentation, registration, and subject-specific atlas estimation[END_REF][START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF], however, the complexity and computation burden for large datasets may limit their use. Regarding intra-subject 4D registration, Lorenzi et al. [START_REF] Lorenzi | 4D registration of serial brain's MR images: a robust measure of changes applied to Alzheimer's disease[END_REF] have proposed 4D non-linear registration via a global 4D deformation optimization scheme in the Demons registration framework.

To compare anatomical differences, 3D population templates have proven their importance for different applications such as mapping (function, structure, vasculature, etc) [START_REF] Thompson | A framework for computational anatomy[END_REF], volume estimation [START_REF] Collins | ANIMAL: Validation and Applications of Non-Linear Registration-Based Segmentation[END_REF] and group comparisons [START_REF] Ashburner | Identifying global anatomical differences: Deformation-based morphometry[END_REF]. While different techniques exist to create unbiased population templates for cross-sectional studies [START_REF] Guimond | Automatic Computation of Average Brain Models[END_REF][START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF], few of these techniques have been developed for the creation of an individual subject 3D template. Reuter et al. created a 3D template for longitudinal analysis by computing the median of the linearly registered subject images [START_REF] Reuter | Within-subject template estimation for unbiased longitudinal image analysis[END_REF].. In this article, a new method is proposed to analyze longitudinal MRI volumes by creating robust 3D individual templates for each subject through non-linear registration and local spatio-temporal regularization of 4D registration. We show that a local spatial constraint over time can have positive global effects to significantly reduce random noise in the measurement of structure volumes such as the lateral ventricles.

Methods

In the following text, we consider all images to be in the same stereotaxic space (i.e., the space of an AD template [START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF]) after going through the same standard preprocessing: non-uniformity correction [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF], brain masking [START_REF] Eskildsen | BEaST: Brain extraction based on nonlocal segmentation technique[END_REF], intensity normalization [START_REF] Nyúl | On standardizing the MR image intensity scale[END_REF] and affine registration [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF]. Similarly to a population template method, individual templates are obtained after averaging the co-registered scans in a hierarchical manner to a subject-specific minimum deformation template [START_REF] Miller | Statistical methods in computational anatomy[END_REF], as described below and illustrated in Fig. 1. The objective of the template creation algorithm is to find the non-linear transformations that minimize the anatomical shape differences between images to create the most representative average of the subject's anatomy. Therefore, to obtain an unbiased individual template, we need to optimize the non-linear transformation with respect to spatio-temporal constraints. While the spatial regularization accounts for smooth and diffeomorphic deformation between each time-point and the template, the temporal regularization enforces continuity in the time domain. In the following sections, the method and the notation for the template creation is inspired from Fonov et al. [START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF], and the nomenclature is presented in Table 1.

Individual template creation

The implementation of the 4D non-linear registration is done in the framework of a 3D non-parametric vector field estimator, subject to certain regularity constraints. The non-linear deformation field is estimated in a hierarchical manner where the algorithm maximizes the local cross-correlation of the image intensity of the source images (i.e., the images ( )

v I i from different time points i, i=[0.
.n]) and the target image (i.e., the current estimate of the template ( ) v Φ ). In this framework, given a set of images ( ) v I i , the individual 3D template must satisfy an intensity constraint (Eq. 1) and a deformation constraint (Eq. 2). The intensity constraint minimizes the intensity difference of the set of images and the current template, ( )

v I Φ such that: ( ) ( ) ( ) ( ) ∑ ∫ = Φ Φ - Φ = Φ n i volume i I dv v I v v i 1 2 * , min arg ) ( ψ (1)
where

( ) v I * Φ
is the optimized subject-specific individual template of subject I, and * ,Φ i ψ is the non-linear transformation that maps the visit i at each voxel location v to current template. The nonlinear registration transformation is also spatially constrained with an elastic body model that minimizes the intensity difference of the paired images (i.e., between template and time point images). The elastic body constraints are justified in such intra-subject registration where large deformations are not expected. The parameters of the elastic regularizer are chosen to ensure that the transformation defined by the vectors is smooth, bijective and invertible.

The deformation constraint removes bias by minimizing the sum of all deformations for all time points:

∑ ∫ = Φ Φ = n i volume i i dv v v 0 2 , , | ) ( | min arg ) ( ψ ψ ψ (2)

Spatio-temporal regularization

The final constraint for the nonlinear registration is performed in the spatio-temporal domain in order to obtain a smooth non-linear deformation over time, since we expect the anatomical changes to happen in a continuous fashion. Therefore, to create a smooth longitudinal deformation field, we perform a local linearization of the deformation fields in the local neighbourhood Ω v and linear regression in time of the zeroth order component ( ) (t T ) and a first order component ( ) (v J ), see Eq.3. The proposed decomposition of the deformation fields is similar to a Taylor series expansion of order 1 in time and in space. The spatio-temporal regularization is effectively performing bi-linear regularization of the local deformation field, see Eq. 4,. Figure 2 provides a simplified diagram of the spatio-temporal regularization process.

( ) ( ) Resulting bi-linear regularization of the deformation field is shown in gray as well as the regularized local Jacobian.

v x t v J t v T t v v x I - + = ∑ Ω ∈ ) , ( ) , ( , * ψ (3)
t v T v T t v T ⋅ + = ) ( ) ( ) , ( ˆ1 0 and t v J v J t v J ⋅ + = ) ( ) ( ) , ( ˆ1 0 (4)

Optimization and convergence

The optimization is performed at 3 hierarchical levels, starting with deformations estimated every 8mm, then every 4mm and finally, every 2mm. At each level, the regularizations are performed consecutively in the order of equations 2 and 3. Also, the template and deformation fields estimated at one hierarchical level are used to initialize the procedure at the next successive hierarchical level. In previous studies, we found that 4 iterations are enough for the convergence of the iterative process at each hierarchical level [START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF].

Experiments

The proposed method, temporal bi-linear regularization with individual template (TWT), is compared to two more classical approaches using the same 3D nonlinear registration framework. The first approach, cross-sectional (CS), is a direct non-linear registration of each time point independently within the common stereotaxic space (AD template) without using any longitudinal information or constraints [START_REF] Collins | ANIMAL: Validation and Applications of Non-Linear Registration-Based Segmentation[END_REF]. The second approach, cross-sectional with individual template (CSWT), creates a subjectspecific template from all time points, but no temporal regularization stage is performed. This will enable us to evaluate the contribution of temporal regularization.

In order to evaluate the stability, regularity, continuity and bias of the proposed approach, we choose to perform ventricular segmentation for each subject at each time point. Ventricular enlargement is often used as a surrogate of brain atrophy, believed to represent the progression of neurodegenerative processes as in MS or AD [START_REF] Nestor | Ventricular enlargement as a possible measure of Alzheimer's disease progression validated using the Alzheimer's disease neuroimaging initiative database[END_REF]. In addition, ventricular segmentation is relatively straightforward due to the high contrast between brain tissue and CSF, and as such should not add unwanted variance to the analysis, as would be the case if more subtle structures such as the hippocampus or entorhinal cortex were used as a test structure.

For each method, the manually segmented lateral ventricles of the AD template are then propagated back to each 3D image for each time-point by concatenation of the required transforms to avoid blurring due to multiple interpolations. Better techniques, more accurate and sensitive to perform ventricular segmentation exist, as described by Apostolova et al. [START_REF] Apostolova | Hippocampal Atrophy and Ventricular Enlargement in Normal Aging, Mild Cognitive Impairment (MCI), and Alzheimer Disease[END_REF], but our goal is to show the longitudinal performance and consistency of our approach.

Data

Two datasets were used to evaluate the proposed algorithm. First, to evaluate stability and potential bias, a scan-rescan database of 20 healthy subjects scanned 4 times within the same week (twice the first session and once during 2 different days). The subject was taken out from the scanner before getting back in for each rescan session.

No ventricular volume change is expected in this case. The T1-weighted MRI images were acquired on a 1.5T SIEMENS MRI scanner with a 3-D spoiled gradient echo (GRE) sequence (TR=22ms, TE=9.2ms, flip=30 o ).

Second, to evaluate the algorithm where change is expected, T1-weighted images from the ADNI study were used. The experiment was performed using all subjects that had 4 visits (0, 6, 12 and 24 months) yielding groups of 177 normal controls (NC), 200 mild cognitive impaired (MCI) and 108 AD patients.

Results

Scan-rescan dataset

Figure 3 shows the temporal distribution of ventricular volume changes between the cross-sectional (CS, and CSWT) and the longitudinal technique (TWT) for the repeated sessions, and where the first session is used as a reference. As expected, the longitudinal regression reduces the variability of volume changes between successive acquisitions. The individual template created with the CSWT approach presents less variability than the segmentations produced with the CS technique. The comparisons of the variances (ANOVA) revealed significant differences between the different sessions only with the CS approach (p-value<0.05).

Longitudinal database from ADNI

Figure 4 shows the percentage of ventricular volume change for each group (NC, MCI and AD) where the first visit is used as a reference (i.e., percent change from baseline is plotted). While the trend of ventricular enlargement is similar between the methods, we can observe the net decrease of variability with the TWT method.

Smaller longitudinal variability should improve the statistical power to detect group differences and thus reduce the number of subjects required for a clinical trial. For example, in order to detect at least a 25% reduction in the annualized rate of ventricular enlargement in the AD group, the TWT approach requires half as many subjects than the standard cross-sectional method (see Table 1).

CS

CSWT TWT Annualized growth rate (%) [START_REF] Thompson | Bias in tensor based morphometry Stat-ROI measures may result in unrealistic power estimates[END_REF].97 (+/-3.65) Table 2. Sample size needed (for treated and placebo groups) to detect a 25% reduction in the annualized rate of ventricular growth in the AD group, with a one-sided test and α = 0.05 at 80% and 90% power using the cross-sectional and the two longitudinal approaches. Note that the TWT method requires roughly half as many subjects when compared to the CS method. Figure 5 shows the TWT-segmented ventricles and the Jacobian determinant of the deformations estimated (from template to time-point, and not from time-point to baseline image) with the 3 methods for an AD patient over the 24-month period. Note that the average Jacobian determinant of the individual visits has been subtracted from the CS Jacobian map to remove the AD template-to-individual deformation and keep only the longitudinal deformation component for comparability with the CSWT and TWT Jacobian maps. Fig. 5. Ventricle segmentations, templates and deformation fields. First row represents the ventricular segmentation performed by the TWT approach. The ventricle increases from 64cc to 70cc in the 24m period. The next row shows the AD template used as a stereotaxic space for CS. The following two rows shows the individual template computed with CSWT and TWT. At each visit, the determinant of the visit-to-template deformation's Jacobian (in blue (larger than template) -red (smaller than template)) is overlaid onto the subject T1W axial image (in grey).

First, we can note that there are large local deformations near the cortex and in some regions of the parenchyma in the CS results, probably due to mismatch with the AD population template (left side) used by the CS method. Furthermore, there are multiple punctuate shrinking and enlarging regions within the ventricles that are not consistent with the notion of gradual ventricular growth. By using a subject-specific template in the CSWT and TWT methods (leftmost image, bottom two rows), there are no large local deformations in the parenchyma or cortex. The CSWT and TWT maps present globally smooth progressive Jacobian changes over time, with much less heterogeneity compared to the CS approach. Surprisingly, the CSWT Jacobians are quite consistent over time. However, there are non-null deformations recovered at the 12 month time point where the individual should be similar to the subject-specific average template.

Discussion and conclusion

In this article, we presented a new approach for the estimation of individual longitudinal changes using unbiased individual subject-specific templates and spatiotemporal regularization. The robust estimation of the deformations is obtained using an unbiased individual template approach, minimizing deformations between subject visits. Meanwhile, a local spatio-temporal regularization is achieved with bi-linear regression of the deformation field and its Jacobian matrices. The regression of the decomposition enables a temporal regularization at a voxel and local neighbourhood level.

Experiments on real data reveal increased stability in estimating individual changes over time compared to standard cross-sectional approaches. The longitudinal regularization of the deformation at a local level reduces the longitudinal noise at the global/structural level, while the hierarchical iterative process produces a robust individual template that allows for better anatomical matching across time in an individual.

In this study, the small number of time points limits the use of more complex longitudinal regression models. In the future other models of temporal regularization regression will be investigated.

An important aspect of longitudinal clinical and research studies is the cost of recruiting and scanning subjects for multiple visits. Adapted longitudinal analysis techniques, such as the one proposed here, will allow for better power to detect differences between groups, and thus will lead to the reduction of the number of subjects required for research and for clinical trails. Conversely, such methods could also shorten the time needed to achieve significant results.

Fig. 1 .

 1 Fig.1. Individual model template algorithm with spatial and temporal regularization. Nonlinear registration is computed pair-wise at each level between the MRI from each visit and the current average template. After convergence, the template is used as a registration target for the next level. Meanwhile, the temporal information of the subject visit allows for spatio-temporal regularization of the deformation fields.

  set of deformation fields, ^ represents the linear regression of each component x d , y d , z d of the deformation vector at voxel v . and first order vectors obtained from the linear regression of the deformation field vectors (blue vectors in Fig 2) and J 0 (v), J 1 (v) are the zeroth and first order matrices obtained by performing regression on the local Jacobian matrix:

Fig. 2 .

 2 Fig. 2. Spatio-temporal regularization. The figure represents a set of 2D longitudinal deformation fields (1 to 4) where the considered vector is red, and the neighbourhood vectors are blue. Resulting bi-linear regularization of the deformation field is shown in gray as well as the regularized local Jacobian.

Fig. 3 .

 3 Fig. 3. Percentage of ventricular volume change for the scan-rescan dataset, measured between the baseline scan and rescans 1, 2 and 3. The use of an individual template with (CSWT and TWT) reduces the bias of the longitudinal registration. The temporal regularization (TWT) results in smaller variability between successive sessions.

Fig. 4 .

 4 Fig. 4. Percentage ventricular change of the ADNI database after 6, 12 and 24 months. The AD group of patient shows a stronger ventricular growth. While TWT estimates smaller ventricular growth compared to CS and CSWT, TWT shows a much smaller variability.

Table 1 . Notation

 1 Deformation fields of subject I at voxel v and at t, time of the visit Jacobian matrix at voxel v of the local neighbourhood Ω

	) ( ) (v v ( ) v i Φ Φ I i ,	Set of images for subject I from different time points i at voxel v Individual template at voxel v
				from baseline
	J	(v	)

ψ Deformation field of visit i to template Φ at voxel v ( ) t v I , ψ

Acknowledgements: Canadian Institutes of Health Research (MOP-111169 & 84360), les Fonds de Recherche Santé Québec and the MS Society of Canada (Bio-Med PhD studentship, 691) and the Alzheimer's Disease Neuroimaging Initiative