Rana Mhanna

Ruchun Zhou

Prof Ming Jin
email: mingjin@tongji.edu.cn

Haiyan Pan

Decheng Wan

Davy-Louis Versace

Nelly Hobeika

Amine Khitous

Olivier Soppera

Patrick Lamielle

Jean-Pierre Malval
email: jean-pierre.malval@uha.fr

Geometry-driven Mass Transport Dynamics within Permeable 3D-Microstructures fabricated by Two-Photon Polymerization

Keywords: Two-Photon Polymerization, Photoinitiators, Micropatterning, Diffusion

The geometric scalability offered by two-photon stereolithography (2PS) has contributed to its rapid widespread into multiple well-established or emerging applications. Taking benefit of this flexible parameterization, we have developed an original fluorescence imaging method enable to map the real-time diffusion of quenching species moving throughout permeable three-dimensional (3D) microstructures with modular geometries (cone, pyramid, square…) fabricated by 2PS. In this strategy, we first methodically characterize the 2P polymerization performances of a Y-shaped triphenylamine photoinitiators series whose initiation efficiency perfectly echoes that measured upon one-photon excitation. More particularly, we show that the photoinitiation mechanism implies a covalent integration of the triphenylamine-based fluorophores into the photopolymer leading to the microfabrication of blue emissive 3D-objects which can be spatially 'switched-off' during the inner diffusion of copper (II) cations used as oxidative quenchers. The modulation of microstructures geometry not only regulates the global symmetry of diffusion profiles but also accounts for a topological control of the quenching propagation dynamics at local scale.

INTRODUCTION.

Two-photon polymerization (2PP) constitutes one of the most promising maskless stereolithography technologies [1][2][3] allowing the production of truly 3D-structures with high spatial resolution. The generated 3D-architectures can display a very high degree of structural details which are tailored specifically using a very flexible computing assisted parameterization which is now routinely proposed in commercially available two-photon 3Dprinting systems 4 . From a fundamental point of view, the exceptional high writing accuracy of 2PP is achieved by tightly focusing of a fs-pulse laser into the photoactive materials. This simple optical configuration promotes an intrinsic spatial confinement of the photopolymerization reaction due to the synergistic combination of a low probable twophoton absorption process scaling quadratically with the excitation intensity and a polymerization threshold effect which can reduce drastically the elementary polymerized volume far below the focus optical dimension 5,6 . With this respect, the typical resolution for near infrared 2PP fabrication is about 100 -200 nm 2,3 which clearly makes this technique an ideal toolbox for the production of intricate 3D-objects from nano-to microscale. In this context, a very large panel of functional materials have been successfully oriented toward 2PP stereolithography including molecularly imprinted polymers 7 , metamaterials 8 , stimuli responsive materials 9 , hydrogels [10][11][12][13][14][15][16] , bactericidal agents 17 , photoclickable dyes 12,18 or magnetic nanoparticles 19 . Such cross-coupling strategies give rise to a new generation of integrated micro/nano devices whose functionalities performance can be modulated and amplified through geometrical and/or topological effects driven by the scale reduction. For instance, Lu et al. 20 have developed a 2PP method generating polymer-based 3D architectures incorporating a large amount of aligned carbon nanotubes whose orientation is both driven by the spatial confinement effect 21 and the laser scanning direction. As a result, it was demonstrated that the electrical conductivity of a linear microstructure fabricated with a parallel laser scanning direction was impressively amplified by three orders of magnitude as compared to that of a similar linear microstructure fabricated through a perpendicular direction. A very versatile strategy to introduce specific functionalities within the final twophoton patterned materials consists in implementing the added value either into the crosslinkable prepolymer or into the two-photon activable initiator. Concerning the first approach, Sun et al. 9 fabricated for instance a biomimetic microlens whose focal distance can be reversibly tuned by changing the pH value. The corresponding 2PP-compatible formulation comprised a bovine serum albumin used as the cross-linkable stimuli responsive protein and the methylene blue as two-photon sensitizer. In the second alternative method, Versace et al. 17 have for instance developed a light-triggered antimicrobial approach based on the 2PP fabrication of 3D polyacrylate microcages with embedded curcumin dyes which led to 95 % bacteria mortality after only 10 min of visible-light exposure. In this highly virulent bactericidal strategy, the 2PP formulation simply consists in a triacrylate monomer associated with the curcumin dye playing a pivotal role as two-photon initiator and as reactive oxygen species photogenerator. In line with this latter methodology, we present herein a multifunctional two-photon initiators series whose 2PP properties and performances are first methodically investigated. We show that the free radical photoinitiation mechanism of these Y-shaped multipolar dyes (see Scheme 1) leads to an efficient covalent integration of the triphenylamine fluorophores into the native two-photon generated macrostructure. By crosscoupling such an unexpected blue emissive property with the permeability of the 3Dmicrostructure, we develop an original fluorescence mapping method to probe at local scale the real-time diffusion of an oxidative quenching reaction whose propagation dynamics can be judiciously modulated with the geometry and topology of the 3D-microstructures.

Scheme 1.

Molecular structures of the chromophores and their electronic level diagrams connected to the excitonic model.

EXPERIMENTAL SECTION.

2.1 Materials. The synthesis and the characterization of all the triphenylamine-based chromophores are described elsewhere 22 . N-methyl diethanolamine (MDEA), copper (II) perchlorate hexahydrate, Rhodamine B, polyethylene glycol diacrylate monomer (PEGDA, Mn 575) and pentaerythritol triacrylate (PETIA) were purchased from Aldrich. All the solvents which are spectroscopic grade were also purchased from Aldrich.

Steady-state absorption and fluorescence spectra.

The absorption and fluorescence spectra were recorded using a Perkin Elmer Lambda 2 spectrometer and a FluoroMax-4 spectrofluorometer.

Cyclic voltammetry and spectroelectrochemistry.

The redox potentials of all triphenylamines were measured by cyclic voltammetry using a Radiometer Voltalab 6 potentiostat. All measurements were performed in N2-saturated acetonitrile with 0.1 M of (nBu)4NPF6 used as supporting electrolyte. The cyclic voltammograms (CVs) were recorded using a three-electrode cell with a platinum disk as working electrode. A platinum wire was used as auxiliary electrode. All cell potentials were measured based on a saturated calomel electrode (SCE) reference which was placed in a separate compartment with a salt bridge containing the supporting electrolyte. The spectroelectrochemical measurements were carried out at the surface of a semi-transparent platinum gauze working electrode using a PEEKbased microcell in three-electrode configuration which is described in reference 23 . Similarly to cyclic voltammetry, the spectroelectrochemical measurements were performed in N2saturated ACN with 0.1 M of (nBu)4NPF6.

Two-photon excitation fluorescence.

The two-photon absorption (2PA) spectra of the dyes in dichloromethane (DCM) were recorded using the two-photon excitation fluorescence (2PEF) method 24 based on femtosecond pulse laser excitation with a mode-locked Ti: Sapphire laser (Coherent, Chameleon Ultra II: pulse duration: ~140 fs; repetition rate: 80

MHz; wavelength range: 680-1040 nm). The technical details relative to the 2PEF set-up are described elsewhere 25 . All 2PA cross-sections (d) were determined using fluorescein in water at pH = 11 as standard 26,27 . The experimental error on the reported d is 15 %.

Free radical one-photon polymerization.

The kinetics for free radical photopolymerization was conducted in situ using a FT-IR Thermo-Nicolet 6700 spectrometer.

A drop of each acrylate-based formulation was sandwiched between two polypropylene films to reduce the oxygen inhibition effect. A 30 µm-thick teflon spacer was also positioned to guarantee a constant thickness of the sample which is then inserted between two BaF2 pellets and irradiated at 365 nm using a Mercury-Xenon Lamp (LC 9588/02A from Hamamatsu) equipped with a 365 nm band pass filter (A9616-07 from Hamamatsu). The photopolymerization rate was determined based on the time-dependent decrease of the vinyl C=C stretching vibration band recorded at 1630 cm -1 .

2.6

Two-Photon Microfabrication. A detailed description of the 3D-stereolithography set-up and the methodology can found in reference 28 . All the microstructures were fabricated upon two-photon polymerization using the same mode-locked Ti: Sapphire oscillator as employed for 2PEF. Each microstructure was fabricated based on a systematic procedure : i)

Deposition of the formulation drop on a 18 x 18 mm glass coverslip ii) Direct laser writing inside the resin iii) Development of 3D-structure by rinsing with ethanol to wash away the unexposed resin.

Polymerization conversion rates of two-photon fabricated microstructures.

The polymerization conversion rates from two-photon patterned microstructures were evaluated by surface-enhanced Raman scattering (SERS) spectroscopy. Multiple 5 x 5 arrays of polymer squares (20 x 20 µm) were two-photon fabricated on SERS-active metal surfaces. The preparation and characterization of these substrates are detailed elsewhere 29 . These surfaces consist in gold nanoparticles deposited on glass coverslips (22 x 22 mm) by a thermal dewetting method from a sputtered gold film. The SERS signals from the polymerized microstructures were recorded using a confocal micro-Raman system (LabRam 300 from Jobin Yvon-Horiba) equipped with a BX40 optical microscope from Olympus. A 532 nm diode-pumped solid-laser (CW) was used as an excitation source. The signals were collected through a 50X, 0.50-NA objective in back-scattering optical configuration by a spectrometer equipped with a thermoelectrically cooled CCD detector. The laser power was fixed at 12.3 mW with an integration time of 20 s for each SERS measurement.

Atomic Force Microscopy (AFM)

. AFM images were recorded using a PicoPlus 5500 atomic force microscope used in ambient atmosphere. The characterization was obtained in resonant mode using a silicon noncontact probe (Arrow Silicon SPM Sensor).

2.9 Fluorescence image microscopy. The time-dependent fluorescence quenching reactions within the 3D-microstructures were monitored using a home-made chamber composed of two PTFE pieces which were designed in order to sandwich the 18 x 18 mm glass coverslips containing the two-photon polymerized microstructures at their upper surface (see Scheme S1A). The PTFE cell was magnetically assembled using a sealing gasket leading to an inner volume of ~ 150 µl. The chamber is then positioned an Olympus IX73 inverted microscope equipped with a 75 W Xe lamp housing. The microstructures immerged in ACN solution were continuously excited from the bottom of the chamber at lexc. = 365 nm through a 40X, 0.65-NA objective using an Olympus U-FUN fluorescence mirror unit composed of a 365 nm band pass filter (BP360-370), a dichroic mirror (DM410) and long pass filter (BA420IF). The epifluorescence images were recorded with a CMOS colour camera. The procedure for fluorescence image monitoring can be described as follows: i) 75 µl of ACN are first introduced into the chamber leading to the swelling equilibrium of the polymer microstructures.

ii) The fluorescence images from the incubated µm-objects are collected upon excitation at 365 nm. To minimize fluorescence bleaching effects, the excitation power is then adjusted using neutral filters placed on the excitation optical path. iii) Fluorescence images collection is then triggered just after the rapid injection of 75 µl of an ACN of Cu 2+ at 8 mM into the chamber. To guarantee a reliable reproducibility of the kinetics, each measurement was repeated with three similar samples (see Scheme S1B).

RESULTS AND DISCUSSION.

Figure 1 shows the one-and two-photon absorption spectra of the chromophores in dichloromethane (DCM). The 1PA spectra of these derivatives have been extensively characterized elsewhere on the basis of excitation anisotropy measurements and TD-DFT calculations 22 . Briefly, each absorption spectrum is dominated by an intensive and structureless band in the 320-430 nm range. The CV of AMK presents a reversible oxidation wave (Eox = 0.95 V vs. SCE) which leads to the generation of a N-centered radical cation whose stability contrasts to that of the triphenylamine (TPA • + ). Indeed, it is well known that TPA • + is unstable due to its rapid dimerization reaction yielding a subsequent tetraphenylbenzidine by-product 36,37 . The rate of this oxidative radical-radical coupling can be drastically reduced by para-substitution effects on the triphenylamine group [38][39][40] . AMK typically exemplifies such a structural change with the introduction of a para p-conjugated substituent (i.e. phenone moiety). The absorption spectrum of AMK • + can be easily observed by spectroelectrochemistry as presented in to the formation of a stable amine-centered radical cation similarly to that observed for AMK whereas the second wave leads to the generation of dication whose second positive charge is presumably located on the external anisole (AMO) or thioanisole (AMS) moiety. It is worth noting that whereas AMO • + (Eox 1 = 1.00 V vs. SCE) is generated at lower potential than 2+ (Eox 2 = 1.28 V vs. SCE) is produced at higher one than AMS 2+ (Eox 2 = 1.23 V vs. SCE). This reverse electrochemical trend suggests the planarization of the AMS biphenyl branch which results in a higher first oxidation potential due to an extent of the p-conjugation but leads to a decrease of Eox 2 consecutive to a lower coulombic repulsion between the two positive charges of AMS 2+ displaying a near planar quinoid-type structure 41,42 . Interestingly, such a planarization is also in line with the bathochromic shift observed for the longest wavelength absorption band of AMS •+ (lMAX ~ 1020 nm) as compared to that of AMO •+ (lMAX ~ 930 nm). A comparable bathochromic shift can be derived when comparing the absorption spectra of the dications (see Figure S1).

AMS • + (Eox 1 = 1.05 V vs. SCE), AMO
Therefore, one clearly shows that the simple change of a methoxy p-substituent by a thiomethyl one significantly affects the biphenyl twisted geometry upon dyes oxidation. Upon light excitation, we have previously demonstrated 22 that the deactivation of the singlet excited state of the dyes mainly proceeds through two competing relaxation processes : i) a radiative deactivation process from an emissive intramolecular charge transfer state ( 1 ICT) ii) an intersystem crossing (ISC) process generating reactive triplet species available for free radical photoinitiation. Interestingly, we showed that the 1 ICT state of AMS was highly polar with a dipole moment twice higher than that of AMO. This distinctive property at the excited state corroborates the consecutive electronic and geometrical changes of AMS during its oxidation as compared to AMO. Within our series, we also showed that polarity of the 1 ICT state indirectly regulates the ISC efficiency with a subsequent impact on the dyes photoreactivity.

This original structure-reactivity relationship gives rise to a clear photoreactivity hierarchy within this triphenylamine series. More precisely we showed that AMO exhibits the highest reactivity when used as type II photoinitiator for free radical photopolymerization. For example, Figure 3 shows the double bonds conversion time-profiles for different triacrylate formulations upon irradiation at 365 nm. Each bicomponent formulation comprises a viscous triacrylate monomer, the pentaerythritol triacrylate (PETIA, see structure in Scheme S2),

which was mixed with the chromophores at isomolar concentration and an aliphatic amine, the N-methyl diethanol amine (MDEA). The first step of the photoinitiation corresponds to a hydrogen abstraction from MDEA to the excited triplet state of the photoinitiator which both generates a a-aminoalkyl radical (MDEA • ) and the chromophore ketyl one. Due to its very high reactivity 43 , MDEA • is mainly responsible for a rapid addition reaction onto the double bonds of the triacrylate monomer. In Table 1 are summarized all the photopolymerization data such as the polymerization rates (Rp/[M]) and the final conversions after 90 s irradiation. (

In this equation, Iabs corresponds to the light absorbed by the sample. kp and kt are the rates of the polymerization propagation and termination steps respectively. These two latter parameters which mainly depend on the nature of the monomer (M) remain constant for our triacrylate formulations. The relative values for Fi are gathered in Table 1. As previously demonstrated 22 , AMO is the highest efficient photoinitiator of the series with an initiation efficiency multiplied by a factor is 11.8 and 3.3 as compared to that of AMK and AMS respectively. Note that these amplification factors measured in highly viscous resin (h25°C ~ 1200 cps) are significantly lower than those measured previously in low viscous diacrylate formulations 22 (h25°C ~ 15 cps). Therefore, it can be derived that the medium viscosity appears Table 1.

One-and two-induced polymerization data relative to the triacrylate resins (Isomolar formulations of photoinitiators: 7 mM / MDEA : 25 mM ).

This hierarchy in the photoreactivity observed for this photoinitiators series can be also extended to the two-photon polymerization. For this purpose, the two-photon initiation efficiencies of the same formulations were determined using a methodology based on the twophoton fabrication of polymer lines upon varying the laser excitation power. Figure S2 typically shows a set of a photopatterned lines written by tightly focusing a fs-pulse laser beam (lexc: 780 nm) into the triacrylate formulation of AMS using a 40X, 0.65-NA objective.

By gradually decreasing the excitation power, it is possible to determine the value of the twophoton polymerization threshold power (Pth) for a given excitation wavelength. This fundamental parameter corresponds to the minimum local absorbed-energy density required to promote the two-photon polymerization (see star pictogram in Figure S2). Interestingly, the value of Pth echoes the minimum concentration of radicals [R0] which ensures the degree of monomer conversion to reach the gelation point at macromolecular scale. Therefore,

[R0] and Pth are connected by the following relationship 2, 45 :

(

Where d.Fi and C respectively correspond to the two-photon initiation efficiency and to the concentration of the photoinitiator in the formulations. It is also assumed that the initiation sequences consecutive to one-or two-photon absorption are similar in such a manner that Fi remains constant for both photoactivation processes. [R0] is also considered as invariant for a given monomer. Therefore, it can be derived from equation 2 that d.Fi is proportional to 1/Pth As shown in Figure 4A, each two-photon polymerization action spectrum matches very well the red edge of its corresponding 2PA spectrum (see Figure 1). This spectral effect clearly indicates that Fi is wavelength-independent upon two-photon excitation. Similarly to onephoton polymerization, AMO exhibits the highest initiation efficiency (Fi ~ 1/d.Pth 2 ) with a relative value at 740 nm which is enhanced by a factor 8.2 and 3.1 as compared to that of and AMS (see Table 1). Interestingly, these amplification factors agree those measured for one photon polymerization as nicely illustrated in Figure 4B where the initiation efficiencies of dyes upon two-photon polymerization at 740 nm (1/d.Pth 2 ) correlate very well their respective initiation efficiencies upon one-photon irradiation at 365 nm (Fi in Table 1).

This correlation clearly confirms that the one-and two-photon polymerization reactions both proceed through an equivalent initiation mechanism.

From a practical point of view, we have previously shown that our triphenylamine derivatives constitute relevant type II photoinitiators for one-photon 3D-printing 22 . In the present case, we can extend their application to the domain of multiphoton 3D-stereolithography. As an illustrating example, we used AMS formulation with PETIA monomer to fabricate µm objects of arbitrary geometries. The microstructures were directly impressed into the triacrylate formulation by 3D translation of the fs-laser focus point. The final object is finally obtained through a development step using ethanol which dissolves away the unexposed monomer. Due to its very efficient cross-linking ability, PETIA monomer leads to a very rigid macromolecular structure which guarantees the mechanical and geometrical integrity of the photopatterned structures. For instance, Figure 5A shows a scanning electron micrograph (SEM) depicting an array of two-photon polymerized cubes spatially organized into a 3 x 3 matrix. On going from the left to right columns of the matrix, each cube of 15 µm side was periodically replicated applying a 1.5 homothetic size reduction factor. The SEM images of a more intricate 3D structure depicting a µm-gear with 26 cogs is also shown in Figure S3. Such a geometrical architecture is typically employed for MEMS applications 46,47 . Two other microstructures are also shown in Figure S4 and µJ. In order to determine the diffusion coefficient for the quenchers, the fluorescence quenching kinetics of each microstructure were analyzed according to a planar diffusion model proposed by Winnik et al. 50,51 . In this formalism which assumes that the timedependent fluorescence quenching is mainly controlled by the progressive z-diffusion of the quenchers along of the microstructure thickness, the emission intensity from the entire structures at time t can be expressed by:

(

where I0 and L denote the initial fluorescence intensity and the microstructure thickness respectively. B corresponds to a parameter connected to the response of the system at the equilibrium state (B = I0/Ieq.-1). Since we consider that all microstructures are initially incubated an ACN solution with a constant concentration of quenchers, the boundary condition appropriate to calculate the integral in equation ( 3) leads to:

In equation ( 4), D corresponds to the global diffusion coefficient of the quencher in the photopolymer. The solid lines in Figure 8B show the best-fitting curves to the experimental data using least squares analysis. According to this methodology, the diffusion coefficient (D)

of the quencher has a calculated value of 2.02 ± 0.04 x 10 -2 µm 2 s -1 for the microstructures generated at 22 µJ whereas D is about 0.79 ± 0.05 x 10 -2 µm 2 s -1 for the square planes fabricated at 51 µJ. Therefore the diffusivity of the quenchers within the two-photon patterned materials can be regulated by the laser writing intensity which indirectly drives the crosslink density of the polymer network. In the same manner, the modulation of microstructures geometry can also drive mass transport dynamics of the quenchers and consecutively the symmetry of their diffusion profiles within these permeable objects. As an illustrating example, Figure 9 shows different time-dependent diffusion profiles within four 3Dmicrostructures which were two-photon polymerized with specific geometries (see insets of programmed by 3D-stereolithography. For instance, the quenching half-time for both pyramids exhibits their lowest value at x = LMAX in line with a very fast quenching dynamics at the thinnest zone of the structure. On going to the center of the pyramids (i.e. x → 0), a significant lengthening of t1/2 is observed as a result of the progressive increase of the thickness (see Figure S10). Interestingly, the presence of a deep hollow in the center of the truncated pyramid can be evidenced not only by the quenching diffusion maps (Figure 9B) but also by the x-profile of t1/2 (Figure 9C). Indeed, one can observe that the quenching half-time reaches its maximum at the middle x-path of the pyramid side then gradually decreases on going to the center. This is clearly in line with a reduction of the pyramid thickness when probing this hollow zone. Therefore, the evolution of t1/2 and thereby the quenching reaction dynamics along the x coordinate clearly accounts for the z-topological of the 3D-object. With a pyramidal architecture, we showed that the quenching propagation occurs from the edge to the center with a t1/2 maximum close to x = 0. As shown in Figure 9C, the diffusion response of two concentric annular structures with different radii results in two t1/2 maxima for x » 0 and » LMAX. As a result, the quenching dynamics alternatively accelerates and decelerates on going from the center to edge of the structure. This dynamics modulation nicely demonstrates a geometry-driven process which can be easily regulated at µm scale by 3D-multiphoton fabrication.

CONCLUSION.

The two-photon polymerization ability of a Y-shaped triphenylamine-based series has been methodically investigated from the dyes electronic properties to the quantification of their two-photon initiation performances. These type II photoinitiators which were initially designed for free radical one-photon polymerization in the UV-Visible range clearly display a transposable reactivity upon two-photon activation in the NIR region. More particularly, we showed that the comparative polymerization performances upon two-photon absorption parallel those measured upon one photon activation which confirms a similar initiation mechanism for both light activation processes. It is noteworthy that this mechanism implies the covalent integration of the fluorescent triphenylamine subunit into the native macromolecules giving rise to a two-photon patternable blue emissive polymer. Taking benefit of these valuable properties, we developed an original methodology using fluorescence imaging microscopy to spatially probe the real-time progress of a fluorescence quenching reaction diffusing throughout permeable 3D-microstructures whose geometry and topology are modulated by two-photon stereolithography. This strategy authorizes not only the local regulation for the reaction dynamics but also the control of the global symmetry associated to the diffusion profiles. This geometry-driven approach can be undoubtedly extended to drug delivery applications with, for instance, the 3D-fabrication of µm-reservoir whose specific geometry can presumably regulate at the local scale the dynamics of drugs release.

Figure 1 .Figure 1 ,

 11 Figure 1.One-(full lines) and two-photon (circles) absorption spectra of the dyes in DCM.

Figure 2 .

 2 Figure 2. A. Cyclic voltamogramms of the compounds in ACN + (nBu)4NPF6 (0.1 M) on platinum electrode at 100 mV s -1 (concentration of chromophores: 10 -3 M). B. Spectroelectrochemical spectra of the triphenlamine derivatives during the progressive generation of their radical cations when applying an oxidation potential at their respective Eox 1 value (solvent: ACN + (nBu)4NPF6 (0.1 M)).

Figure 2B .

 2B Figure 2B. The radical cation in ACN exhibits an intensive band in NIR region with a maximum at 760 nm. This band is clearly redshifted with respect to that of the tri-ptolylamine radical cation 40 (DlMAX ~ 90 nm) which confirms the extension of the charge delocalization along the biphenyl branch. As shown in Figure 2A, the electrochemical behaviours of AMO and AMS are similar. Both derivatives present reversible CVs which comprise two successive one-electron oxidation waves. The first oxidation step corresponds

Figure 3 .

 3 Figure 3.

  which also drives the output of the balancing mechanism between ICT and ISC relaxation processes. a Iabs ≈ I0.el.[PI], b Rp / [M] : maximum of the first derivative of the conversion rate time profile, c Relative value which comprises the constant factor kt/kp 2 (see equation 1)

Figure 4 .

 4 Figure 4. A. 2PA action polymerization spectra of the triacrylate formulations. B. Linear correlation between one-and two-photon polymerization efficiencies of the formulations. (l1PA = 365 nm and l2PA = 740 nm).

Figure 5 .

 5 Figure 5. A. SEM images of representative 3D cubes structures fabricated by 2P induced polymerization. (PETIA with AMS / MDEA : 7 mM / 25 mM, lex : 780 nm, P = 5 mW, texp = 2 ms). B-C Transmitted light microscope image and corresponding fluorescence image of a 2D logo using the same acrylate-based formulation and the same microfabrication conditions (Scale bars: 10 µm).

Figure S5 .

 S5 The first one presents a hollow icosahedral architecture positioned on a circular polymerized pedestal. Such a 3D microstructure which is more sophisticated than the previous one is also more fragile since its lateral edges are formed by suspended thin polymerized lines (1 µm diameter) which can be deformed due to surface tension during the final washing process. The second hollow structure in FigureS5exhibits an interesting waffle-like topology composed by multiple 20-µm-high vertical walls with a length of 50 µm and a very thin thickness of 450 nm which leads to a high-aspect-ratio of 44. Finally, a last example consisting in a 2D logo is depicted in Figures5B-C. Interestingly, this logo is highly fluorescent when excited at 365 nm. This blue emission presumably stems from the excitation of AMS by-products whose integrity of the p-conjugated triphenylamine core is preserved during the photoinitiation step. Moreover, this triphenylamine subunit should be covalently connected to polymer network due to the addition of AMS ketyl radical to acrylate double bonds. To demonstrate the occurrence of this parallel reaction, FigureS6depicts the progressive change of AMS fluorescence spectrum during the photopolymerization of the triacrylate formulation at 365 nm. The photopolymerization reaction induces both a total collapse of AMS fluorescence ICT band located at lMAX = 480 nm and the concomitant rise of a new emission band in the high energy region (lMAX = 390 nm). The significant band blueshift observed between the fluorescence spectra of AMS and its by-product corroborates the reduction of AMS ketone functions consecutive to the addition of the ketyl radical to the acrylate C=C bonds. It should be noted that an alternative radical-radical coupling process is less likely here because the resulted by-product would not be covalently linked within the polymer. To demonstrate such an effective grafting, we have prepared a new acrylate-based formulation associating the same bi-component system (i.e. AMS and MDEA) mixed with a poly(ethyleneglycol)diacrylate monomer (PEGDA, see Scheme S2). Contrary to PETIA resin, photopolymerization with PEGDA formulation leads to a permeable polymer which allows mass transport of embedded species when incubating the polymer into an imbibing solvent. According to this methodology, it is therefore possible to evaluate whether the byproducts are grafted within the photopolymer or can be leached out from the macrostructure.

Figure

  Figure S7depicts the evolution of the fluorescence signal recorded from a polymerized dot during its incubation in acetonitrile (ACN). This object was fabricated by one photon polymerization at 365 nm with PEGDA resin and corresponds to an impressed half focus point using a 0.65-NA objective. As shown in FigureS7, the polymer dot radius increases by a factor ~ 1.31 in ACN which confirms the solvent uptake during the polymer swelling. In parallel, the fluorescence signal from the polymer dot as well as its fluorescence image remains perfectly invariant during 20 min incubation in ACN. Therefore, it can be derived that all the fluorophores are immobilized within the photopolymer presumably due to their covalent linkage within the matrix. Hence, the photopolymerization of AMS formulations leads to a photopatternable polymer that becomes emissive due to a concomitant covalent integration of the fluorophore into the native macrostructure. Such an unexpected emission property can be oriented for the real-time mapping of a well-chosen quencher diffusing within the permeable microstructure whose dimensions and geometries can be easily tuned by 3D-

Figure 6 .

 6 Figure 6. A. Optical transmitted image of the photopolymerized cone incubated in ACN. B. SEM images of the cone (Scale bar: 45 µm) C. Time-dependent fluorescence images of the polymer cone in ACN upon progressive diffusion of Cu 2+ within the polymer network. D. Evolution of the fluorescence profile during the diffusion quenching process. Inset: Plot of the radius at fluorescence quenching half-time (R1/2) as function of the diffusion time of Cu 2+ within the polymer cone.

Figure

  Figure 6D presents various normalized fluorescence profiles for distinctive diffusion times.The propagation rate of the quenching front can be derived from the time-dependent change

Figure 7 .

 7 Figure 7. A. SEM micrographs of four polymer square planes fabricated with distinctive thicknesses. (Scale bar: 20 µm) B. Corresponding time-dependent fluorescence quenching images of the excited squares during their incubation in the same ACN solution with 4 mM Cu 2+ . C. Log-log plots of the squares thicknesses vs. their corresponding quenching half-times.

Figure

  Figure 7B displays fluorescence images for each square at different diffusion times.According to the symmetry of the microstructures, a z-diffusion process (i.e. perpendicular to the plane) is clearly privileged. For instance, the time-dependent decrease of images lateral size is ~ 3-times lower as compared to that observed along z-coordinate. This dominant z-diffusion process also follows a Fickian regime as nicely illustrated by the square-root relationship connecting planes thicknesses with their corresponding quenching half-times (see Figure7C). It is worth noting that the stability and the morphology of the two-photon polymerized squares remain globally invariant after the quenching reaction. This is typically illustrated in FigureS9which presents the AFM images of a dried square plane whose fluorescence emission was quenched upon incubation of the structure in an ACN solution with 4 mM Cu 2+ . Note that the shrinkage of the polymer after drying the microstructure gives rise to a surface roughness with an arithmetic average value of about 120 nm. Interestingly, the diffusivity of the quenchers within these microstructures can be modulated with the writing laser intensity. Indeed, such a photonic parameter can be adjusted in order to regulate the conversion rate of the acrylates C=C bonds and thereby the corresponding crosslink density of the two-photon patterned materials. With this respect, a set of two distinctive 20 x

Figure 8 .

 8 Figure 8.A. SERS spectra in the 1520-1690 cm -1 range for two 20 x 20 µm square planes which were two-photon polymerized using two distinctive laser excitation energies. All SERS spectra were normalized by the intensity of 1730 cm -1 (nC=O mode) B.Normalized fluorescence quenching kinetics of the corresponding excited squares during their incubation in an ACN solution with 4 mM Cu 2+ . The solid lines correspond to the best-fitting curves using equation 3.

Figure 9C )

 9C Figure 9C) : i) A square plane of 20 µm side and 1.5 µm thick. ii) A square base pyramid of 20 µm side and 20 µm height. iii) A similar pyramid was fabricated but with a truncation at h = 10 µm. The upper side of this truncated pyramid comprises a conical-shaped hollow with 5 µm radius and 7 µm deepness. iv) A 1 µm-thick polymer disk associated with two concentric annular structures of 3 µm height and with radii of 3 and 10 µm respectively. The quenching reaction progress within each microstructure has been monitored by subtracting the fluorescence image at a given diffusion time to the initial fluorescence one (i.e. I0 (x,y) -It (x,y)). According to this methodology, the real-time propagation front of diffusing quenchers can be mapped over the entire structure as depicted in Figure 9B.

ACKNOWLEDGEMENTS.

This work was supported by the National Natural Science Foundation of China (51873154, 51573139).

ASSOCIATED CONTENT.

CORRESPONDING AUTHORS.

* M. Jin (e-mail: mingjin@tongji.edu.cn), J.P. Malval (e-mail: jean-pierre.malval@uha.fr)

AUTHOR CONTRIBUTIONS.

NOTES.

The authors declare no competing financial interest.


Vincent Hakim
email: vincent.hakim@ens.fr

Simon Pigeon
email: simon.pigeon@lkb.upmc.fr

Amandine Aftalion
email: amandine.aftalion@ehess.fr

Metamorphosis of the Landau transition in the flow of a resonantly-driven bistable polariton fluid

The onset of dissipation in a flowing superfluid has been the subject of numerous investigations. Motivated by recent experiments, we theoretically analyze the analogous phenomenon for a onedimensional "quantum fluid of light" which is resonantly driven, and exhibits bistability. The flow is found to abruptly change multiple times when the fluid velocity or the potential strength are increased. In contrast to the classical Landau transition to a time-dependent flow, the transitions take place between stationary states and involve the fluid bistability in an essential way.

It is well-known since the classical analysis of Landau that superfluids exhibit a critical velocity at which dissipation sets in. The discovery of Bose-Einstein condensation in cold atomic vapors has allowed detailed studies of superfluidity [START_REF] Pitaevskii | Bose-Einstein condensation and superfluidity[END_REF]. The Landau transition and defect emission in the flow past an obstacle have been experimentally observed with the creation of vortices [START_REF] Raman | [END_REF], or gray solitons in a one-dimensional setting [3]. It is welldescribed theoretically in the framework of the Gross-Pitaevskii equation [4][5][6][7][8]. Beyond cold atomic vapors, exciton-polariton fluids are attracting significant attention [9], as solid-state devices with a higher condensation temperature allowed by the exciton-polariton very low mass. Bose-Einstein condensation has been achieved in these "quantum fluids of light" [9][10][11] and superfluidity has been studied numerically [12] and experimentally [13,14], including the nucleation of vortices at a critical flow velocity in the wake of an obstacle [14]. Polaritons are composite bosons that result from the strong coupling between the excitonic resonance of a semiconductor quantum well and a microcavity electromagnetic field [9]. In early experiments, the driving field needed to create polaritons was used either transiently [13] or in a spatially localized way [14] to avoid interfering with the superfluid behavior. The short polariton lifetime then restricts the experiment duration or limits the observations to a local region around the pumping spot. In order to bypass these limitations, it is useful to introduce a support field away from the strong pumping spot used to create the polaritons [15]. This new coherently driven regime has started to be investigated theoretically [11,[15][16][17][18] and experimentally [19][20][21][22]. The quasi-resonant drive tends to lock the phase of the condensate but when the support field is not too strong, it nonetheless allows the formation of vortices [20] and dark solitons [21,22]. However, the dynamical properties of this driven-dissipative condensate and their dependence on the fluid bistable character [23] remain to be better understood [15,[20][21][22].

Here, we analyze the flow past an obstacle of a resonantly forced condensate in the bistable regime [11,15,17,[20][21][22]. We focus on the one-dimensional case which is the easiest to theoretically analyze. We find, as for conventional superfluids, that a sudden flow change can occur in the wake of the obstacle when the fluid velocity is increased or, when the obstacle strength is increased at fixed velocity. However, the characteristics of this transition and of the resulting flow are found to be quite different from the Landau transition of a superfluid and depend in an the essential way on the fluid bistability.

We consider the fluid described by the following generalized Gross-Pitaevskii equation (GGPE)

i ∂ t ψ = -2 2m ∂ 2 x ψ+ V (x) -∆ -i γ 2 + g|ψ| 2 ψ+F e ikpx .
(1) In the context of exciton-polariton microcavity physics, Eq. ( 1) provides an effective description of a driven lower polariton field [9,11], with the polariton-polariton repulsive interaction accounted by the constant g > 0. Additional terms as compared to the usual GPE arise from the coherent drive and dissipation [9,11]. Quasi-resonant pumping is characterized by its amplitude F , its momentum k p , produced by a slight tilt of the driving laser beam with respect to the cavity plan and by the detuning ∆ between the frequency of the driving field and the bottom of the lower polariton band. Dissipation is described by the rate γ > 0 arising from the polariton finite lifetime. The potential V (x) models a localized obstacle. It is our main aim to characterize its effect on the fluid flow described by Eq. ( 1). It is worth noting that our results are also relevant for nonlinear optics where Eq. ( 1) is known as the Lugiato-Lefever equation [24] and describes wave evolution in a cavity filled with a nonlinear medium (see e.g. [25] and ref. therein).

The explicit x-dependence in Eq. ( 1) can be eliminated by defining

ψ = γ/2g φ(x) exp(ik p x) (2)
The function φ then obeys the equation

i∂ τ φ = - 1 2 ∂ yy φ -ik 0 ∂ y φ -δ(y) + i -|φ| 2 φ + f (3)
where we have introduced the dimensionless variables y = x mγ/2 , τ = tγ/2, and constants,

k 0 = k p 2 /mγ, f = F √ g(2/ γ) 3/2
, and defined the function

δ(y) = δ 0 - 2 γ V (y), with δ 0 = 2 γ ∆ - ( k p ) 2 2m (4)
Before considering the effect of a localized obstacle, we briefly recall some properties of the fluid described by Eq. (3). When V (y) = 0 and δ(y) = δ 0 , Eq. ( 3) has constant solutions in space and time with a homogeneous density ρ = |φ| 2 which can readily be seen to simply obey,

B(ρ) := [(ρ -δ 0 ) 2 + 1]ρ = f 2 (5)
Two cases can be distinguished. When δ 0 < 0, the function B(ρ), defined in Eq. ( 5), is obviously monotonically increasing from 0 to ∞ with the density. As a consequence, the density ρ is an increasing function of the forcing amplitude f . When δ 0 > 0, B(ρ) can be non-monotonic with multiple homogeneous solutions for a given forcing. A simple analysis of Eq. (5) shows that this actually happens when δ 0 > √ 3. An example of this S-like dependency of the density with the driving field is plotted in Fig. 1a-insert and leads to bistability in a window of intermediate forcing strengths i.e. for B(ρ + ) < f 2 < B(ρ -) with ρ ± = 1 3 2δ 0 ± δ 2 0 -3 . This bistability for sufficiently strong blue detuning is well-known in nonlinear optics [9] and has been demonstrated for polaritons in microcavities [23]. While in the low-density (LD) regime self-interactions are unimportant, the high-density (HD) regime with strong selfinteraction allows the observation of superfluidity [26] in a polariton fluid.

Having recalled the basic features of the homogeneous state, we proceed and describe our simulations of Eq. ( 3) with a localized repulsive (V (y) > 0) gaussian potential

V (y) = γ 2 u m exp[-(y/σ) 2 ] ( 6
)
We focus on the bistable parameter regime with δ 0 > √ 3 and the forcing f in the appropriate intermediate interval (see Fig. 1a). The fluid density is set on the HD state upstream of the obstacle. In an experimental setting, this would result from a strong driving field in a far upstream local region, as proposed in [15] and experimentally realized in e.g. [21]. For a low potential amplitude, the flow is steady, with the density decreasing as expected in the region of the repulsive potential, and returning smoothly to the HD state in the wake of the obstacle, as shown in Fig. 1a. An increase in the potential amplitude u m produces a transition in the flow, as shown in Fig. 1b, reminiscent of the usual onset of dissipation in a flowing superfluid. However, the character of the transition appears to be very different from the classical Landau transition. Instead of the repeated emission of gray solitons in one dimension [3,5], the flow restabilizes after a transient. It results again in a stationary flow but with the fluid density in the LD state downstream of the obstacle, as shown Fig. 1b. In other words, for the drivendissipative GGPE, the critical velocity signals a metamorphosis of the steady flow instead of the appearance of a time-dependent flow. That the flow density downstream of the obstacle lies in the LD state provides a first hint that fluid bistability is playing a significant role in the observed transition. This is further suggested by the fact that, as in previous works [15,21], a transition is not observed when the driving amplitude f is strong enough for the fluid to be in the monostable HD regime i.e. when f > f + in Fig. 1a. Simulations performed for potential amplitudes very close to the transition (Fig. 1c,d) allow one to better see how the metamorphosis of the steady flow pattern is taking place. When the amplitude of the potential is close to, but below, the critical potential amplitude, the obstacle is followed by a fluid region close to the intermediate density (ID) unstable state (Fig. 1c). This region terminates by a front that joins the ID state to the more downtream HD state. As the critical potential amplitude is varied and the critical amplitude is approached, this front stands farther and farther downstream from the obstacle, with an increasing region of the fluid downstream of the obstacle in the unstable ID state. For potential with an amplitude slightly greater than the critical one, the complementary process is observed, as shown in Fig. 1d. As for subcritical potentials, the obstacle is followed by a fluid region in the unstable ID state but which terminates by a front joining it to the stable LD state. With increasing potential amplitude, this front stands closer to the obstacle. It reaches the obstacle and disappears, as soon as the potential amplitude departs by a small amount from the critical one. An important conclusion from these observations is that the critical solution is such that the fluid downward wake exactly stands at the unstable ID state.

In order to more fully understand this transition and the role of bistability in the stationary flow metamorphosis, we consider the parameter regime suitable for theoretical analysis, provided by an obstacle that varies on a long length scale, σ 1 (Eq. ( 6)). For a slowly varying obstacle, when the flow is accordingly slowly varying, the derivative terms in Eq. ( 3) can be treated perturbatively. At the lowest order, they can be entirely neglected and the "adiabatic" solution, φ a (y) = ρ a (y) exp[iθ a (y)], readily obtained. The fluid density ρ a , is linked to the potential amplitude by Eq. ( 5) with δ 0 simply replaced by δ(y) (Eq. ( 4)), which takes into account the influence of the potential on the detuning. In this adiabatic approximation, solving the quadratic Eq. ( 5) for δ(y) provides the implicit relation between the fluid density and the potential,

δ(y) = ρ a (y) ± f 2 ρ a (y) -1, 0 ≤ ρ a ≤ f 2 . ( 7
)
As for homogeneous solutions, the solution phase is sim-ply given as a function of the density

θ a (y) = arctan 1 ρ a (y) -δ(y) (8)
The relation (7) between the density and δ at fixed forcing amplitude f , is plotted in Fig. 2a. It is equivalent but more convenient for our purpose than the insert of Fig. 1a, which gives the density as a function of f for fixed "detuning" δ. A simple calculation shows that Eq. ( 7) determines the density as a unique function of δ(y) when f < f c = (4/3) 3/4 1.2408 while for f > f c , there is a range of δ values with multiple possible densities. In other words, bistability occurs for a range of δ values when f > f c , as illustrated in Fig. 2a.

Let us now consider, a fluid injection in the HD state, when the forcing is sufficiently strong for bistability to occur (i.e. f > f c ). As the potential varies with the position y, δ(y) follows it according to Eq. ( 4). The density, as given by Eq. ( 7), moves along the HD branch in Fig. 2a. The adabatic solution is already a close approximation of the flow obtained by numerically solving Eq. ( 3), for the gaussian potential of Eq. ( 6) even with a rather large amplitude (u m = 6) when σ = 4 (Fig. 2b).

We first briefly describe the case of an attractive potential (u m < 0). Eq. ( 7) predicts that the fluid density goes up the high density branch as V (y) becomes more negative. The flow should undergo a transition if u m is large enough for the top of the high density branch to be reached, since the branch cannot be followed beyond its top. This transition is indeed confirmed in numerical simulations of Eq. (3) (Fig. S1). As for repulsive potentials, when |u m | is larger than the critical one, the flow is stationary and, above the transition the fluid density in the wake of the obstacle is in the LD state.

How a transition can happen for a repulsive potential (u m > 0) is less obvious. The density of the adiabatic solution (Eq. ( 7)) follows the high density branch toward low density before increasing again in the wake of the obstacle, as shown in Fig. 2a. This appears a smooth process for all potential strengths u m . It is not clear why this would result in a transition of the flow profile at a critical amplitude u m and what this critical amplitude would be. In particular, no hint is provided by computing linear (i.e. Bogoliubov) excitations [18,19] around the local homogeneous fluid density. However, one can note that in the adiabatic approximation, all derivatives are absent and, as a consequence, the fluid velocity plays no role. This suggests to go beyond the adiabatic approximation and treat perturbatively the derivatives terms. Writing ρ(y) = ρ a (y) + ρ 1 (y) + • • • , θ(y) = θ a (y) + θ 1 (y) + • • • , the first corrections to the adiabatic solution of Eq. ( 7), ( 8) are obtained after a short calculation as,

ρ 1 (y) = - 2k 0 B (ρ a ) dρ a dy , θ 1 (y) = k 0 B (ρ a ) d[δ(y) -2θ a ] dy (9)
where the function B(ρ) is defined in Eq. ( 5). These corrections are shown in Fig. 2b and, as expected, they result in a closer agreement between the analytic approximations and the numerical profiles. More interestingly, the corrected density profile in the (ρ, b) diagram provides a clue to the origin of the instability (Fig. 2a). One observes that the correction (9) produces a departure of the profile from the high density branch towards the middle unstable branch when the potential returns to 0, in the close downward wake of the obstacle. Eq. ( 9) shows that this non-adiabatic effect grows with k 0 and, it also grows with the localized potential amplitude u m . One can therefore guess, that, for sufficiently large u m or k 0 , this leads the flow profile loop in Fig. 2a, to reach the unstable density branch in the (ρ, b) diagram and lead to an instability. While suggestive, this argument is obviously not rigorous since the perturbative correction ( 9) cannot be trusted when it provides a significant correction to the lowest order result.

In order to obtain a full reduced nonlinear description, a further asymptotic limit is needed, beyond that of a slowly varying potential (i.e. σ → ∞). A simple mathematical one is obtained by increasing the flow velocity k 0 at the same time as the length scale of the potential is varied, i.e. taking the limit, σ → ∞, k 0 → ∞ with a fixed ratio κ = k 0 /σ. Determining the steady solution of Eq. (3) reduces in this limit to solving the simple system,

κ∂ z ρ = -2ρ -2f √ ρ sin(θ) (10)
κ∂ z θ = [δ(z) -ρ] -f cos(θ)/ √ ρ (11)
where z = y/σ. Eq. ( 10), (11) simply give back for κ = 0 the adiabatic solution (7), (8) and, perturbatively for small κ, the correction (9). But, in the asymptotic limit considered, κ can now take any value. The reduced system (10), (11) has only first-order derivatives in z. Its phase-plane can be easily determined (Fig. S2). Eq. ( 10), (11) have only one spatially-divergent modes from the unstable ID at z = +∞. A simple shooting method determines the amplitude of the potential for which this divergence vanishes and the solution tends at z = +∞ toward the ID branch, as illustrated in Fig. 3a. It brings however a surprise : for given driving parameters, multiple transitions are found by increasing the localized potential amplitude. The fluid density in the wake of the obstacle is in the HD state for low potential amplitudes. At a first critical amplitude, it jumps to the LD state, as described above. When the potential amplitude is further increased, a second transition is found, at which the fluid density jumps back to the HD state. Further transitions are found for still higher values of u m . The loci of these transitions are plotted in the (k 0 , u m ) parameter plane in Fig. 3b. These multiple solutions and the asymptotics of the u * m (k 0 ) branches can be understood by considering Eq. ( 10), (11) for large κ. In this limit, the evolution of ρ and θ with z is slow, except for the potential term δ(z) that evolves on a z-scale of order 1. The density of a solution of Eq. ( 10), (11), that starts in the HD state at z = -∞, does not significantly change when it encounters the localized potential while its phase rotates by an angle ∆θ,

∆θ = - 2 γκ +∞ ∞ dz V (z) = - √ π u m /κ (12)
where the second equality holds for the Gaussian potential (6). For the solution to end up at z = +∞ in the ID state, the phase turn ∆θ has to bring it precisely, on one of the two entering separatrices of the ID unstable state, as shown in Fig. S2. Namely, ∆θ should be equal to ∆θ 1,2 -2nπ, n = 0, 1, • • • with θ 1 -1.844, θ 2 -4.464 for the parameter values of Fig. 3c. The angle values with n ≥ 1 correspond to the solution phase making full rotations before reaching one of the two separatrices. The double series of critical potential amplitudes follows from Eq. ( 12),

u m = -κ (θ 1,2 -2nπ)/ √ π, n = 0, 1, 2, • • • (13)
Eq. ( 13) shows that asymptotically u m depends linearly on κ = k 0 /σ. The slopes given by Eq. ( 13) for the three lowest branches are displayed in Fig. 3b together with the numerically obtained solutions. The two lowest transition branches merge at κ 1.30. The 4th and higher branches cross and recombine at intermediate κ values, producing the bifurcation diagram shown in Fig. 3b . The bifurcation diagram of Fig. 3b holds for the reduced asymptotic problem (Eq. (10,11)). We have checked that away from this asymptotic limit, similar diagrams are obtained for the GGPE with our reference parameters (Fig. S3). For potential range σ = 1 or 2, three transitions are observed upon increase of the local potential strength u m for large enough k 0 (e.g. k 0 = 6.2 for σ = 2) and no transitions are observed for small enough k 0 [START_REF] Eq | for the GGPE in the large k0 and large um limit for potentials of finite range[END_REF] .

In conclusion, we have found that the drive provided by a constant support field deeply changes the Landau transition. The stationary flow profile undergoes a metamorphosis instead of becoming time-dependent. Moreover, for given flow and pumping parameters, successive transitions exist at a discrete number of potential amplitudes. We have shown that these phenomena can be described and understood analytically in suitable asymptotic regimes. These results certainly suggest a careful reexamination of the Landau transition in higher dimensions. We also hope that they will motivate experimental studies of the phenomenon. Finally, we cannot help but wonder, whether the extended switches of the fluid density wake induced by a localized obstacle could provide useful applications in all-optical technology and devices [27]. 3) is plotted together with the adiabatic solution alone (dotted magenta; Eq. ( 7),( 8)) or with the first-order corrections (dashed magenta; Eq. ( 9)). The parameters are σ = 4, um = 6, k0 = 2.75, f = 3.2, δ0 = 6.2.) 13) are shown (dashed). In the different parameter regions, it is indicated whether the flow in the wake of the obstacle tends toward the high (H) or the low (L) density states with the numbers corresponding to the additional 2π dephasing (Eq. ( 13)) of the large κ solutions in their downstream wake, as compared to the H0 and L0 ones. At finite κ the different regions of the same type with different numbers (e.g. H0 and H1) are delimited by nodal lines of solutions (not shown) with vanishing density at a point allowing the required phase jump. The zigzagging transition line at κ 1 has only been computed up to u * m = 24. The other parameters are f = 3.2, δ0 = 6.2. m,1 (dashed green line) and um = 21.0 < u ⇤ m,2 (dashed-dotted green line). It is also shown a circle of radius equal to the modulus of the HD point (solid black line) centered at the origin (solid black circle) as well as its intersections (red diamond) with the entering separatrices of the ID point. The phase di↵erence ✓1 and ✓2 between these intersection points and the HD point are indicated. They provide the asymptotic slopes of the di↵erent transition branches (see Eq.( 13) and the main text). The other parameters are f = 3.2, 0 = 6.2.

FIG. 1 :

 1 FIG.1: Numerical simulations of Eq. (3). (a) Insert : fluid density ρ vs. forcing f as described by Eq. (5). In the parameter regime considered, there are three homogeneous steady states, stable HD (solid blue circle) and LD states (solid green square) as well as an unstable ID one (solid red diamond). (a)& (b) The fluid is injected in the HD state. The curves represent the fluid density at different times and have been shifted with time to highlight the stationarity of the flow. (a) for a small potential amplitude um = 2, the flow is steady and is in the HD state in the wake of the obstacle at y 0. (b) For a larger um = 6, the flow is still stationary, but is in the LD in the wake of the obstacle. (c)&(d)Detail of the transition. (c) um = 5.50 (d) um = 5.501 (see the main text for a discussion). Other parameters are δ0 = 6.2, f = 3.2, k0 = 2.75 which corresponds to the typical experimental values ∆ = 0.5M eV, γ = 0.1M eV, 2 /m = 1M eV µm 2 , kp = 0, 616µm -1 , √ gF = 0.036(M eV ) 3/2 . The potential range is σ = 1 corresponding to 4.5µm, the chosen unit length.

FIG. 2 :

 2 FIG. 2: Adiabatic approximation and corrections. (a) Plot in the (ρ, b) diagram. (b) Plot of the density and the phase vs. y. The stationary numerical solution (solid blue) of Eq. (3) is plotted together with the adiabatic solution alone (dotted magenta; Eq. (7),(8)) or with the first-order corrections (dashed magenta; Eq. (9)). The parameters are σ = 4, um = 6, k0 = 2.75, f = 3.2, δ0 = 6.2.)

FIG. 3 :

 3 FIG.3: Reduced asymptotic description. (a) Solutions of Eq.(10,11) for κ = 2.75. In the obstacle wake, the flow tends towards the HD state for um = 4.9 (dashed blue) and toward the LD state for um = 5.0 (dashed-dotted green). The critical flow corresponds to u * m 4.93 (solid dark) and tends to the ID state. (b) Diagram of the transition lines in the (κ, um) plane (solid black). The asymptotic slopes for large κ of the three lowest transition lines given by Eq. (13) are shown (dashed). In the different parameter regions, it is indicated whether the flow in the wake of the obstacle tends toward the high (H) or the low (L) density states with the numbers corresponding to the additional 2π dephasing (Eq. (13)) of the large κ solutions in their downstream wake, as compared to the H0 and L0 ones. At finite κ the different regions of the same type with different numbers (e.g. H0 and H1) are delimited by nodal lines of solutions (not shown) with vanishing density at a point allowing the required phase jump. The zigzagging transition line at κ 1 has only been computed up to u * m = 24. The other parameters are f = 3.2, δ0 = 6.2.

1 FIG

 1 FIG.S2: Phase plane analysis of the reduced asymptotic problem. (a) Phase plane of the homogeneous problem(Eq. (10,11)) (i.e.without the localized potential) showing the three fixed points, HD (blue solid disk), ID (red solid disk), LD (green solid disk) together with the two entering separatrices (red solid lines) of the ID fixed point and the outgoing ones, ending on the HD point (blue solid line) or the LD point (green solid line). Several trajectories are shown ending either on the HD point (dashed blue line) of the LD state (dashed green line). (b) Same diagram as (a) but showing two critical trajectories of the reduced problem for  = 6.0 with the localized gaussian potential of amplitude u ⇤ m,1 ' 7.328 (dashed red line) and u ⇤ m,2 ' 21.34 (dashed-dotted red line), corresponding to the first two transitions when um is increased from 0. The critical trajectories start in the HD fixed point and they end on the ID fixed point. For each of these transition points, two other trajectories are shown which end either at the HD point, with um = 6.5 < u ⇤ m,1 (dashed blue line) and um = 21.8 > u ⇤ m,2 (dashed-dotted blue line), or at the LD point, with um = 8.0 > u ⇤ m,1 (dashed green line) and um = 21.0 < u ⇤ m,2 (dashed-dotted green line). It is also shown a circle of radius equal to the modulus of the HD point (solid black line) centered at the origin (solid black circle) as well as its intersections (red diamond) with the entering separatrices of the ID point. The phase di↵erence ✓1 and ✓2 between these intersection points and the HD point are indicated. They provide the asymptotic slopes of the di↵erent transition branches (see Eq.(13) and the main text). The other parameters are f = 3.2, 0 = 6.2.

Acknowlegments. We would like to thank A. Bramati and E. Giacobino for instructive discussions.

Supplementary figures

Metamorphosis of the Landau transition in the flow of a resonantly-driven bistable polariton fluid Vincent Hakim, Simon Pigeon and Amandine Aftalion Supplementary figures transition with um = 14.16 (LD state) and um = 14.17 (HD state). Note that the 2nd transition is inverted as compared to the 1st, namely the fluid density jumps back to the HD state when um is increased. While the density at y 1 in the HD states are the same in (b1) and (b2), the phase is lower by 2⇡ in (c2) as compared to (c1) in agreement with the asymptotic analysis. The phases of the LD state at z 1 also di↵er by 2⇡ between (c2) and (c1) while for larger k0 they would the same as in the asymptotic analysis. This additional phase jump comes from the crossing of a (nodal) line of LD solutions, for which the density vanishes at a point, when increasing um in the LD region. (b3) &(c3) 3rd transition with um = 14.84 (HD state) and um = 14.85 (LD state). The transition is analogous to the first one, but with the phases shifted by 2⇡ in the far downstream wake of the obstacle.


M Cassier

T Degiovanni

S Guenneau

F Guevara Vasquez

Active exterior cloaking and mimicking for the heat equation

. We point out the second method can be applied to waves in dispersive media.

I. INTRODUCTION

The temperature u(x, t) of a homogeneous isotropic body satisfies the parabolic heat equation

∂u ∂t = σ∆u + h, for t > 0, (1)
where σ is the thermal diffusivity (a scalar constant), h is a source term and ∆ is the Laplace operator in 2D.

If one makes a change of variables in [START_REF] Cassier | Active Thermal Cloaking and Mimicking[END_REF], σ∆u becomes ∇ • σ ∇ , with σ an anisotropic heterogeneous parameter and ∇ the gradient in the transformed coordinate system. Cloaking of heat in the transient regime was originally introduced in this way that led to an anisotropic heterogeneous diffusivity [START_REF] Guenneau | Transformation thermodynamics: cloaking and concentrating heat flux[END_REF]. An experimental validation of this concept with a metamaterial cloak engineered by the group of Martin Wegener at the Karlsruhe Institute for Technology can be found in [START_REF] Schittny | Experiments on Transformation Thermodynamics: Molding the Flow of Heat[END_REF]. Changing the shape of the cloak and reducing its thickness leads to extreme material parameters (except for diamond shape cloaks). Many research groups have contributed to the design of thermal metamaterials, notably to make them more amenable to fabrication and this field of transformation thermodynamics is now quite mature. We recommend the research monograph for a state of the art [START_REF] Huang | Theoretical Thermotics Transformation Thermotics and Extended Theories for Thermal Metamaterials[END_REF]. However, managing thermal fields does not necessarily require the use of metamaterials, as we shall now see.

II. EXTERNAL CLOAKING FOR THE HEAT EQUATION

We wish to find sources for the heat equation (1) in a homogeneous medium so that they generate a cloaking field u c (x, t) that cancels out a probing field u i (x, t) inside a region Ω without perturbing the total field outside Ω. By linearity of the heat equation the total field u c (x, t) + u i (x, t) vanishes inside Ω allowing us to hide e.g. a homogeneous Dirichlet inclusion, while it is identical to u i (x, t) outside of Ω, preventing an observer outside of Ω from detecting the inclusion solely based on thermal measurements. Of course, our approach is not limited to Dirichlet inclusions, it applies to any linear passive object characterized by a boundary condition (Dirichlet, Neumann, Robin, ...) or by transmission conditions (e.g. penetrable inclusions) or even to some active objects. It is based on the Green identities to show that a distribution of monopoles and dipoles supported on ∂Ω can provide the desired cloaking field [START_REF] Cassier | Active Thermal Cloaking and Mimicking[END_REF]. In practise, these monopoles and dipoles can be realized using heat pumps such as Peltier devices, see Figure 1. This cloaking method inspired by [START_REF] Miller | On perfect cloaking[END_REF] in the context of propagative phenomena comes with two main drawbacks: (a) the incident field u i (x, t) needs to be known ahead of time and (b) the sources needed to generate the cloaking field completely surround the object, as shown in Figure 1. However, this method provides also an efficient way to create thermal illusions, as shown in Figure 2, without dealing with complex metamaterial parameters as in passive mimicking methods [START_REF] Alwakil | Inverse heat mimicking of given objects[END_REF].

III. ACTIVE EXTERIOR CLOAKING FOR THE HELMHOLTZ EQUATION WITH COMPLEX WAVENUMBERS

To work around drawback (b), we notice that the Fourier-Laplace transform of a solution to the heat equation must solve the Helmholtz equation with a complex wavenumber. To be more precise, let u(x, t) be a solution to the heat equation [START_REF] Cassier | Active Thermal Cloaking and Mimicking[END_REF]. Under mild-assumptions, the Fourier-Laplace transform u(x, ω) of u(x, t) is given by

u(x, ω) = ∞ 0 dt exp[iωt]u(x, t), (2)
and satisfies the Helmholtz equation with a wavenumber k given by the relation k 2 = iω/σ:

∆u(x, ω) + (iω/σ)u(x; ω) = -h(x, ω) , (3)
with ω the angular frequency.

Once in the frequency domain we can achieve active exterior cloaking for the Helmholtz equation by generalizing the approach in [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF] from real positive wavenumbers to complex wavenumbers and then take the inverse Fourier-Laplace transform to recover time domain sources for the heat equation that satisfy our requirements. The approach in [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF] consists in using the Green identities for the Helmholtz equation to determine a density of monopoles and dipoles on ∂Ω that achieves cloaking. Then the boundary ∂Ω is split into N dev segments ∂Ω 1 , . . . , ∂Ω N dev and the monopole and dipole density in segment ∂Ω j is "moved" to a new location x j . This can be done using the Graf addition formula and what we obtain are N dev multipolar sources supported at the new locations x j . In the case k > 0, the series defining the multipolar sources are proven to converge [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF] to the desired fields outside of the union of the gray disks in figure 3.

A numerical result of thermal cloaking is given in figure 4. We proved [START_REF] Cassier | Active exterior cloaking for the 2D Helmholtz equation with complex wavenumbers and application to thermal cloaking[END_REF] that the same result holds for complex wavenumbers (excluding the negative real axis) by establishing that the tail of the series can be bounded by the tail of a geometric series with ratio depending only on the point of interest and its position relative to the new source locations x k and the boundary segments ∂Ω j . Such bounds on the truncation error allow us to predict quantitatively the cloaking effect. We also apply the maximum principle to equation (3) to guarantee that the error from truncating the multipolar sources series on a region reaches its maximum on the boundary of the region. IV. CONCLUSION

Fig. 1 :

 1 Fig.1: Example of an arrangement of Peltier elements that could be used to cloak objects (e.g. a kite) inside a twodimensional region Ω, illustrated here by a disk within a heat conducting plate. Each Peltier device is represented by two adjacent red and blue boxes, where the heat flux it can create is oriented in the direction normal to their interface. We have represented both Peltier elements that transport heat within the plate (across the boundary ∂Ω of Ω) and also between the exterior and the plate. This figure comes from the article[START_REF] Cassier | Active Thermal Cloaking and Mimicking[END_REF].

Fig. 2 :

 2 Fig. 2: Numerical example of a kite mimicking a flower object, both with homogeneous Dirichlet boundary conditions. A snapshot the scattered field from the original object at a given time ts appears in (a). In (c), the scattered field from the original object is made to appear as the scattered field from the object in (b) from the perspective of thermal measurements outside of the cloaking region Ω. This figure comes from the article [1].

Fig. 3 :

 3 Fig. 3: Example of surface ∂Ω split into four sections ∂Ω 1 , . . . ∂Ω 4 . where monopoles and dipoles on ∂Ω j are moved to a new location x j . The disks in gray denote the region where the cloaking field u c diverges.

Fig. 4 :

 4 Fig. 4: Temperature distribution at a given time t resulting from (a) a point source in the presence of a kite object with homogeneous Dirichlet boundary condition and (b) same as in (a), but with the cloaking devices activated. The linear color scale ranges from blue to red, with zero being represented in white.

We have proposed a new path to thermal cloaking, based on active sources. Our approach, when performed in the frequency domain, relies upon a Helmholtz equation with a complex wavenumber, also valid for wave problems in dispersive media. Our proposal is a source of inspiration for other groups [9].