Vaibhav Sharma

Elise Azar

Andre P Schroder

Carlos M Marques

Antonio Stocco
email: stocco@unistra.fr

Active colloids orbiting giant vesicles

published or not. The documents may come

Introduction

Interaction of viruses, bacteria and other colloids with biological membranes plays a determinant role in many important processes such as microbial infections, drug delivery and toxicity from nanomaterials. 1 These processes are related to endocytosis, translocation and adsorption of living and artificial colloids in membranes. In biological cells, many passive and active (fuel-consuming) mechanisms 2 occur in the membrane during the recruitment of a particle inside the cell. Such mechanisms are often coupled together and are far from being understood today as shown by the COVID-19 pandemic. 3 Beside membranes, motion of active and passive particles in many cases occurs close to solid walls, fluid interfaces and different kinds of obstacles. 4 These obstacles not only act as a passive confinement but may also interact with the moving particle. As a function of the gap distance between the particle and the interacting interface, active particle motion can be strongly modified with respect to the motion inside an unbounded fluid. Hydrodynamic and long range surface forces together with specific interactions of the two interacting objects may promote new particle dynamics such as transport of cargo, 5 orbital motion, 6 reflection, 7 constrained 2D 8,9 and 1D directional motion. 10,11 repulsion between the particle and the obstacle, which ensures that the particle cannot penetrate inside the obstacle. 14 The flow far from the particle can be characterized without considering local effects due to the swimming mechanism details. In this case, the flow field is given by 𝑢(𝒙) = 𝛼𝑆 𝐷 (𝒙 -𝒙 𝟎 ; ê) + 𝑂(|𝒙 -𝒙 𝟎 | 3 ), where 𝒙 𝟎 and ê are the position and the orientation of the microswimmer respectively.

𝑆 𝐷 is the symmetric force dipole, and 𝛼 is a coefficient. 𝛼 > 0 describes the flows generated by a "pusher" microswimmers: the flow is directed outward on the front and the back of the particle and inward on the side, which leads to an attraction between the swimmer and an obstacle located on the side. A threshold obstacle size, above which a persistent orbital motion occurs was always predicted by these models.

Particle diffusion was also taken into account to describe the particle escape from the obstacle orbit. 14,15 Even in the absence of a spherical obstacle, the emergence of circular motion of self-propelled particles has been reported in various scenarios in the recent literature. 16,17,18,19,20,21,22,23,24 Here, we present experimental results on the dynamics of active Janus particles of spherical shape interacting with giant unilamellar vesicles (GUVs) made of phospholipid bilayers, which are an established biomimetic system for biological membranes. Both particles and GUVs are sedimented on the bottom substrate. To investigate the interaction between an active particle and a GUV, we have focussed our attention on the role of the particle activity in generating a self-propelled motion and on passive mechanisms such as the particle Brownian diffusion.

Material and Methods

Active colloids.

We used bare silica microspheres of radius RP = 1.98 ± 0.05 µm and fluorescent Melamine Resin (MF) particles with RP = 1.15 ± 0.06 µm (microParticle GmbH, Berlin, Germany). To fabricate Janus colloids, first we prepare a monolayer of silica/MF beads on a thoroughly cleaned silica wafer by drop casting (particle concentration is 0.2% by volume). After completely drying, a thin platinum layer is deposited by a metal sputtering technique (Edwards Auto 306 Evaporator). 25 The resulting particles are half coated with a thin platinum layer of 6 ± 1 nm thickness. This thickness was measured by light reflectivity (Multiskop, Optrel, Germany) on the corresponding planar surfaces. [START_REF] Azzam | Ellipsometry and Polarized Light[END_REF] Janus colloids are released from the wafer using 10 µL of glucose solution (0.1 M) along with gentle agitation using the pipet tip. As reported previously, 16 the fabricated particles may present a variability in shape of the Janus boundaries and coating (see Supporting Information, Figure S1).

2.2. Vesicle formation. Giant unilamellar vesicles, GUVs, are prepared using a PVA (Polyvinyl alcohol) gel assisted formation [START_REF] Weinberger | [END_REF] which allows to fabricate vesicles with a high yield. PVA gel is prepared by dissolving PVA in PBS (Phosphate-Buffered Saline,10 g/L) solution at 5% concentration. PVA gel is spread on a PTFE (Polytetrafluoroethylene) plate and dried for 30 minutes at 80°C after. 10 µL of a 99:1 (molar) mixture of POPC (1-Palmitoyl-2-oleoylphosphatidylcholine) and POPC-NBD (POPC fluorescently labelled with Nitrobenzoxadiazole) lipids in chloroform (1 g/L) are spread on the PVA gel and vacuum dried for 15 minutes. The lipid system is then hydrated with 200 µL of sucrose (0.15 M) and allowed to grow for 2-3 hours. The vesicle suspension is then collected and sedimented in 1 mL of glucose (0.15 M) solution. The slight density mismatch between the sucrose solution inside the vesicle and the sucrose/glucose solution in the aqueous medium allows the vesicles to sediment at the bottom of the cell without strongly deforming the vesicle. A viscosity = 1.1 mPa.s of the sucrose/glucose solution was measured by passive microrheology.

Microscopy and Tracking.

Bright-field and fluorescence microscopy was used to observe Janus colloids and giant vesicles in aqueous media. 10 µL of a Janus particle dispersion and GUV solution along with 120 µL of glucose is filled in the sample cell. In all experiments with active colloids, a final concentration of 2% H2O2 (which is the catalytic fuel of the Platinum Janus colloids), was used. Enough time (few minutes) is allowed to pass before observing the active motion of the colloids. At this hydrogen peroxide concentration no significant change of the GUV size distribution was observed (see Supporting Information, Figure S2) and we did not observe any signature of hydrogen peroxide gradients in our experiments. The optical setup is a Nikon Eclipse TE2000 microscope (×60 objective) equipped with CMOS camera (Orca Flash 4.0, Hamamatsu). Videos were recorded in a 10-100 frame per second (fps) range. Particle tracking was done by the open source software Blender v 2.8. For the tracking of the vesicle centre of mass, ImageJ and Python routines were also used.

The obtained trajectories [t (s), x (µm), y (µm)] provide the particle/vesicle centre of mass position in time in the laboratory coordinate system.

Results and Discussion

Orbital motion of Janus active colloids around isolated giant vesicles

For diluted enough particle and vesicle concentrations, we were able to observe the interaction between single active Janus colloids and single GUVs. The colloids are sedimented on the substrate, along which both passive and self-propelled motion occur. GUVs are also sedimented on the substrate without any strong adhesion between the solid wall and the membrane. A 2% H2O2 aqueous solution was always used to trigger the self-propulsion of the silica-platinum Janus colloids. 28,29 In most of our experiments, whenever the trajectory of the active particle reached the vesicle boundary, the active particle initiated a striking orbital motion around the GUV. Figure . 1A shows a typical trajectory of an active colloid orbiting around a GUV observed (see Supporting Information, Video S1). For this particular example the GUV has a larger but comparable size to the active particle of radius RP ≈ 2 µm. For times t < 1 s, the GUV is still not perturbed by the active particle. For 1 < t < 2 s, the silica face of the Janus particle pushes on the GUV membrane, which starts to translate a distance L with approximately constant speed VGUV = 0.7 µm/s. For t > 2 s, the active particle performs a persistent orbital motion around the GUV. Times spent by the particles in the GUV orbits vary from few seconds up to a thousand seconds. If enough time is allowed to pass, the active colloid eventually leaves the GUV orbit and initiates again an active "free" motion before encountering another GUV. Even if the Janus colloid seems to adhere to the GUV, with the Janus boundary oriented radially with respect to the GUV centre (Figure . 1A), no membrane deformation (partial membrane wrapping) 1 or adhesion was observed by fluorescence microscopy. Nanoscale membrane deformations below the resolution of our methods may however occur. By measuring the vertical positions of the GUV and particle equators (by detecting their respective focal planes using the z-stage of the microscope), we confirmed that both colloid and vesicle sediment at the bottom of the container, and the active motion of the Janus colloids occurs in the wedge defined by the planar substrate and the GUV (see Figure 1B). Orbital trajectories were observed for both of the particle sizes that we used (MF, RP ≈ 1 µm and silica, RP ≈ 2 µm) and in the whole range of GUV sizes investigated here. We observed both rotation directions with a comparable statistics between clockwise (CW) and counter clockwise (CCW) orbital motion directions.

For RP ≈ 1 µm active colloids, in some cases, we also observed switching events between the two rotation directions during an ongoing orbital motion: an active colloid rotating CW (CCW) is able to remain close a GUV and change to CCW (CW) direction, without escaping the GUV, see Figure 2 (two CCW-CW and one CW-CCW switching events). This switching may result from both in plane and out of plane Brownian rotational motion of the colloid.

For RP ≈ 2 µm active colloids, instead, we did not detect switching events between the two possible orbital directions, see Figure 2.

From the circular path of particle trajectories, one could evaluate an orbital radius RO (see Figure 1B and 2) and estimate the size of the GUV. Assuming spherical objects in close contact resting both on the planar substrate, RGUV ≈ RO 2 /(4RP). For some trajectories we also measured directly the GUV size (in fluorescence or bright field microscopy) by changing the focal plane from the particle equator to the GUV. These GUV size measurements agree with the estimated value from the orbital radius. Preliminary Reflection Interference Contrast Microscopy (RICM) experiments show in some cases the signature of a vesicle flattening close to the substrate due to the gravity, which is however very small compared to the vesicle size.

Trajectories shown in Figure 2 demonstrate that orbital motion occurs in a wide range of GUV and active particles sizes, covering vesicle to particle radii ratios RGUV/ RP ≈ 1.5 to 20. Note that some trajectories in Figure 2 are clearly affected by some drift velocity due to convective flows. However, this drift velocity is significantly lower than the particle speed and it can be easily decoupled from the particle motion in our analysis. We observe that active colloids display persistent orbital motion around GUVs, that can be viewed as soft spheres with fluctuating surfaces, a feature also seen for active rod particles orbiting around isolated solid obstacles; and for spherical active Janus colloids orbiting in a crystal made of large solid particles. 6,12 In the bulk, far-field hydrodynamics is able to explain such persistent orbital motion for active particles behaving as a "pusher" squirmers. This squirmer category pushes the liquid in the front and in the rear, and entrain the liquid on the sides (here the Janus boundary), which leads to an effective hydrodynamic attraction with an obstacle (here the GUV) located in close proximitiy. 6,14 However, our experimental results for sedimented active colloids and GUVs show also a clear difference with respect to the bulk far-field hydrodynamic predictions for orbital motion around either spherical solid or liquid obstacles. 6,14 In fact, in our experimental system we did not observe any critical radius of the GUVs for the onset of the orbital motion (see Supporting Information, Figure S3). Hydrodynamic models predict that in the bulk, for obstacle size lower than a critical value RC (≈ 10-15 RP), the active particle shows a scattering behaviour; whilst for obstacle radii larger than RC, a persistent orbital motion of the active particle around the obstacle occurs.

Analysing about two hundred trajectories showing the interaction between an active colloid and a GUV, we evaluate the time spent by the particle around the GUV, which we call the orbital time (torb).

In Figure 3, we plot the normalized counts of orbital time torb (or trapping time 14 ) for both particle sizes.

It varies from few seconds to several minutes. An exponential function exp(torb/t * ) can be used to approximately describe the torb distribution. 6 From the fits, we find characteristic times t * = 23 ± 1 s for RP ≈ 1 µm and t * = 46 ± 5 s for RP ≈ 2 µm (see inset in Figure 3). It is important to remark that t * for RP ≈ 2 µm particles is approximately twice that of RP ≈ 1 µm particles. Interestingly, this proportionality correlates thus with the bulk translational friction coefficient of spheres 6RP. It is thus tempting at this stage to speculate on the possible reasons for a particle to escape the orbit where it has been captured.

The strength of the attraction of a pusher to a wall is known to depend on the velocity of the particle. 10 A significant slowdown of the orbital velocity could then allow for escaping the orbit. However, all the escaping events that we have observed were not associated with a stalling or a significant reduction of the particle velocity. Another possible explanation for orbital escape could reside on the spontaneous stochastic rotation of the particle into an angular position where its propulsion would drive it away from the membrane. 14 Indeed, as we will further study below, the orbital motion of the Janus pusher proceeds with a propulsion direction close to the tangent of the circular orbit, a strong deviation from this orientation could potentially drive the particle away. However, we have in a few instances observed clockwise movement spontaneously changing into counterclockwise, and conversely. During such changes of orbital direction, the particles are favorably oriented to move away from the wall but stay trapped instead. Hence, the escape dynamics may be governed by the Brownian translational diffusion of the particle as predicted by Spagnolie et al.; 14 or by the membrane-particle repulsive interaction observed in the absence of H2O2. Indeed, our passive colloids do not tend to adhere to GUVs in thermal Brownian motion, even when colloids were very close to the GUV membranes.

Active colloid-GUV interaction

In order to elucidate the strength and the nature of the interaction between an active colloid and a giant unilamellar vesicle, we performed a detailed analysis of the trajectories and of images taken during the orbital motion and compared them to data describing active colloids far away from the GUVs. From the active particle tracking, we calculated the discrete velocity vector 𝑣̅ = (𝑣 𝑥 , 𝑣 𝑥 ) where 𝑣 𝑥 (𝑡) =

<𝑣̅ 𝑖 (t+t) • 𝑣̅ 𝑖 (t) > = 4Dtr(t) + V 2 cos(ωt) fr , (1)
where the first term 4Dtr(t) is related to the two-dimensional translational Brownian diffusion, Dtr;

and the second term V 2 cos(ωt) fr, describes the active propulsion. If the active motion occurs in an unbound fluid (i.e. "free" motion), fr = exp(Drt ), and the bulk Brownian rotational diffusion, Dr, completely randomizes the persistent active motion at long time scales. Equation (1) can be used for interpreting free trajectories even in the presence of experimental convection. Indeed, for t > t1,

VACF of the velocity 𝑣̅ can be written as the sum of the active propulsion term in equation ( 1) and a constant drift velocity Vd (see Figure 2): 9

<𝑣̅ (t+t) • 𝑣̅ (t)> = V 2 cos(ωt) fr + Vd 2 . ( 2
)
The orbital motion of the particle around a GUV is not in general described by equation (1). The interactions of the particle with the vesicle introduces a number of other force sources. Beyond the particle propulsion (translation and rotation) and Brownian noise leading to the form in equation ( 1):

dynamical confinement forces and the resulting confined Brownian rotation and translation motions, vesicle shape fluctuations, all conjugate into a complex evolution of the particle orbital trajectory. It is however simple to show that for a perfectly circular orbital motion at constant orbital frequency ω, the autocorrelation function would have the form of equation ( 1) with fr =1, thus providing a convenient manner to directly extract orbital velocities from the particle trajectory data alone.

In Figure 4A, we show a typical trajectory of an active colloid orbiting around a GUV for about seven times (Orbit1). Then, the particle escapes from the first orbit and perform "free" active motion before encountering a second GUV (Orbit2). VACF corresponding to this trajectory divided in three ranges (Orbit1, Free and Orbit2) are shown in Figure 4B.

In the orbital and "free" motion, <v could be also fitted by equation ( 2) with a zero angular velocity (no orbital motion), which also provides fr = exp(Drt ) and leads to the rotational diffusion in the "free" motion. Fitting the data by equation 2, we were able to evaluate, V and ; and measure a negligible drift velocity (< 0.15 µm/s). We evaluate  in the range 0.25 -0.76 rad/s for RP ≈ 1 µm and in between 0.15 and 0.5 rad/s for RP ≈ 2 µm and verified that consistently V =  Rt for all orbital motions. Moreover, rotational diffusions in the free motion Dr,free = 0.10 ± 0.03 rad 2 /s for RP ≈ 1 µm and Dr,free = 0.03 ± 0.01 rad 2 /s for RP ≈ 2 µm were obtained by the fits of equation 2, which are comparable to the calculated bulk values Dr,bulk = kBT/(8RP 3 ) = 0.11 rad 2 /s (RP ≈ 1 µm) and Dr,bulk = 0.02 rad 2 /s (RP ≈ 2 µm).

In Figure 5, we plot the particle speed in the orbit as a function of the speed when the active particle is far from the GUV (for each particle whose trajectory allows a robust measurement of <v(t+t) v(t)> = V 2 in the "free" and orbital motion, see Figure 4B). Even if the H2O2 fuel concentration is fixed at 2%, a range of velocities is observed due to different observation times (from few minutes to some hours after the addition of H2O2, which is consumed in time) and an intrinsic variability in the catalytic coating due to the Janus colloid fabrication. 16 Most of the data lies on the line Vorb=Vfree, in agreement with results obtained for active rods orbiting around solid spherical obstacles, and pointing to a speed independent from the active particle-obstacle distance. 6 However, for the velocities above 3 µm/s, several points show Vorb< Vfree. The equality of velocities in the free and orbital motions can be rationalized within two opposite scenarios. Obviously, velocity equality would hold if no additional dissipation was involved as the particle is captured into the orbital motion. However, this is unlikely due to the proximity to the membrane wall and to the strong distortion of the flow environment in the confinement wedge. Vorb=Vfree could also result from an increase of the propulsion mechanism due to a wall effect, 10 which would be compensated by an increase of the translational drag. 32 Whether this perfect compensation is a coincidence or a consequence of the self-propelled motion confined wedges remains to be elucidated. In order to gain a deeper insight into the rotational dynamics of an active colloid bound to an orbital trajectory, we performed image analysis of both the particle and the GUV around which it orbits. For some cases of RP ≈ 2 µm particles, we were able to detect both the orientation of the particle Janus boundary and the centre of mass of the GUV (see Supporting Information: Video S2 and Video S3). In this way, we could find the particle location angle  and orientation angle  in the GUV frame, see Figure 6. Hence, to evaluate the Brownian rotational diffusion of the Janus particle in the orbit we calculated the mean squared angular displacement < 2 >, which is plotted as a function of the lag time in Figure 6B.

For very short lag times, < 2 > shows a linear behaviour; whilst a saturation could be observed already before t = 1 s, which is related to a confinement effect as also seen in Figure 6A. Such saturation is the signature of confined rotational diffusion due to the fluid flow generated by the active particle, which tends to align the particle Janus boundary towards the GUV centre. 14 For confined Brownian dynamics, < 2 > can be written as: 10

< Δ(𝜑 -𝜃) 2 >= 𝐷 r,orb Γ [1 -exp(-2Γ Δ𝑡)] , ( 3
)
where  is an effective elastic constant for the angular velocity:  =  which tends to restore the particle angular position  (see inset in Figure 6B), with the particle Janus boundary orientated radially to the centre of the GUV. from the fit shown in Figure 7. Hence, < 2 > data confirm a slowing down of the rotational diffusion with respect to the bulk rotational diffusion, which can be associated to a friction due to the active colloid-GUV interaction.

Force and torque transfers can indeed occur when an active particle interacts with a GUV. Already in Figure 1, we showed that a GUV can translate a distance L with a speed VGUV = 0.7 µm/s when an active colloid comes in close contact. An effective force FGUV on the GUV can be calculated as FGUV = tr VGUV ≈ 6RGUV VGUV ≈ 46 fN, where tr is a simple estimate for the GUV translational friction coefficient. These few tens of femto-Newtons were measured several times in our experiments during the impact of active colloids on GUVs.

An additional experimental result is now presented to describe how an active colloid is able to impart also a torque on a GUV. Figure 7 shows a complementary observation for an active colloid (RP ≈ 2 µm), which is not translating but only rotating in a CCW direction around an axis parallel to the normal to the substrate and passing through its centre of mass. Presumably, the colloid is partially stuck onto the solid substrate, but keeps its rotational activity. A GUV (whose size is approximately twice the particle size) revolves around the active Janus colloid also in a CCW direction. The GUV position appears to follow the orientation of the Janus boundary, with the centre of mass of the GUV slightly off with respect to the radial direction defined by the Janus boundary. Angular velocities of the GUV in the orbital motion and Janus colloid rotation are approximately the same,  0.155 ± 0.005 rad/s: a passive GUV is orbiting around an active colloid with a tangential velocity VGUV = (RGUV + RP) ≈ 0.9 µm/s! An effective tangential force acting on the centre of mass of the GUV can be calculated as FGUV = tr,GUV VGUV ≈ 72 fN. Some defects present on the GUV (see Figure 7) also allow monitoring the rotation of the GUV around its own axis. The GUV rotates in a CW direction (opposite to the orbital motion) with GUV = 0.415 rad/s, and a corresponding torque M = r,GUV GUV  8RGUV 3 = 41 fN.µm, where r is the GUV rotational friction coefficient. The simultaneous force and torque acting on the GUV, generated by the local fluid flow of the active colloid, could be explained considering that the centre of mass of the GUV is slightly off with respect to the Janus boundary direction. An effective lever arm l = M/ FGUV = 0.57 µm could thus describe the offset between the GUVs centre of mass and the location of the effective force acting on the GUV (see Figure 7, 40 s).

To summarize, our experiments (Figure 1

Conclusions

We have experimentally investigated at the single particle level, the interaction between an isolated active colloidal particle and a giant unilamellar vesicle. Persistent orbital motion was observed in the whole range of GUV-particle size ratios investigated. Hence, no critical size for the onset of orbital motion was detected, which is in net contrast to far field hydrodynamic models for solid or fluid spherical obstacles in the bulk. In our experimental system, the tensionless and fluctuating nature of lipid bilayers leads to hydrodynamic boundary conditions that may be different to the ones considered so far. 14,15 We also clearly show that a transfer of forces and torques occurs during the active colloid-GUV interaction, which results in a persistent translation of the GUV during the particle impact, a slowing down of the particle rotational diffusion with respect to the bulk, and a persistent GUV orbital motion around a rotating but not translating active particle. Still, far-field hydrodynamics can be used to explain not only the orbital motion due to a fluid flow generated by the active particle and pointing towards the particle Janus boundary, but also the role of the translational Brownian diffusion in the escape dynamics. 14,15 Finally, our results provide new perspectives on the interaction between active particles and soft compartments, which go beyond the effect of guiding 11 and scattering 33 as observed for the interaction of a bacterium with a passive colloid. Active particles are indeed able to impart forces and torques on giant lipid vesicles, as it was very recently reported for active particles encapsulated inside a vesicle. 34 These results motivate future investigations in partial and complete particle engulfment regimes to study directional transport of GUV by active colloids and activity trigged particle endocytosis.

Figure 1 .

 1 Figure 1. (A) Centre of mass trajectory of an active colloid of radius RP ≈ 2 µm interacting with a GUV. Time lapse images of the interaction between an active particle and a GUV are shown in the insets. (B) Sketch of the experimental geometry: a Janus colloid orbiting a GUV along the substrate. z is the difference between the focal plane of the GUV and the particle equator. RO is the radius of the orbital trajectory (see also next Figure).

Figure 2 .

 2 Figure 2. Centre of mass trajectories of active Janus colloids (RP ≈ 1 and 2 µm) around GUVs of different sizes. The time range of the observations is 30 -788 seconds. Clockwise (CW) and counterclockwise (CCW) directions are also reported close to the trajectories.

Figure 3 .

 3 Figure 3. Normalized counts of orbital times of RP ≈ 1 and 2 µm particles rotating around GUVs. The red and black curve are exponential functions exp(torb/t * ). A log-linear representation of the data is shown in the inset.

  (t+t) v(t)> as a function of the lag time is approximately constant, which shows the absence of three-dimensional motion and negligible out of plane rotational dynamics of the active Janus colloid. A decay of the autocorrelation function of the active velocity modulus v = √𝑣 𝑥 2 + 𝑣 𝑦 2 could in fact indicate a non-zero z-component of the velocity, 𝑣 𝑧 . 31 <𝑣̅ (t+t)• 𝑣̅ (t)> in Orbit1 shows a clear oscillation due to the angular velocity  In Orbit1, damping of <𝑣̅ (t+t)• 𝑣̅ (t)> does not prevent the evaluation of the angular velocity. During the "free" motion, a decay of <𝑣̅ (t+t)• 𝑣̅ (t)>

Figure 4 .

 4 Figure 4. (A) Trajectory of a RP ≈ 2 µm active colloid orbiting around two GUVs and performing a "free" motion in between. (B) Velocity autocorrelation functions as a function of lag time Δt for the trajectory shown in (A).

1 Figure 5 .

 15 Figure 5. Orbital speed of RP ≈ 1 µm and 2 µm active colloids as a function of the speed far from the orbit.

Figure 6 .

 6 Figure 6. (A)   as a function of the experimental time for a RP ≈ 2 µm active colloid orbiting a GUV. (B) Mean squared angular displacement < ) 2 > as a function of the lag time. Solid line is a fit of equation 3 to the data. In the inset, a sketch of a Janus particle (P) orbiting around a GUV together with the polar location angle  and the orientation angle  in the GUV frame are displayed.

  , 6 and 7) clearly show the possible translation and rotation of GUVs under the effects of the fluid flows generated by the active particle. Note that our findings are in contrast to previous observations reporting no torque transfer and no rotation of solid spherical particle during the orbital motion of active rods.6

Figure 7 .

 7 Figure 7. Time lapse images of an RP ≈ 2 µm active colloid rotating in a CCW direction (yellow arrow) but not translating and interacting with a GUV, which is both revolving around the Janus colloid and rotating in a CW direction (red arrow). Arrow arc length shows the rotation around its own axis between two frames (time is given on the top of the image).

Acknowledgements

We thank Marc Schmutz for his help with the metal sputtering technique and the Electronic Microscopy platform of the Institut Charles Sadron. E. A. thanks the Emergence@INC of CNRS for financial support. V. S. acknowledges the Doctoral School ED 182 of University of Strasbourg for financial support.

A self-propelled Janus colloid performs a persistent orbital motion around a giant unilamellar vesicle, even when the vesicle size is comparable to the particle size.


Dr Virgile Ayzac

Dr Mohammed Dirany

Dr Matthieu Raynal

Dr Benjamin Isare

Dr Laurent Bouteiller
email: laurent.bouteiller@sorbonne-universite.fr@polymers_su

Energetics of Competing Chiral Supramolecular Polymers

Keywords: bis-urea, helix reversal, net helicity amplification, pathway complexity, supramolecular polymer

Chirality can have unexpected consequences including on properties other than spectroscopic. We show herein that a racemic mixture of bis-urea stereoisomers forms thermodynamically stable supramolecular polymers that result in a more viscous solution than for the pure stereoisomer. The origin of this macroscopic property was probed by characterizing the structure and stability of the assemblies. Both racemic and non-racemic bis-urea stereoisomers form two competing helical supramolecular polymers in solution: a double and a single helical structure at low and high temperature, respectively. The transition temperature between these assemblies, as probed by spectroscopic and calorimetric analyses, is strongly influenced by the composition (by up to 70°C). A simple model that accounts for the thermodynamics of this system, indicates that the stereochemical defects (chiral mismatches and helix reversals) affect much more the stability of single helices. Therefore, the heterochiral double helical structure predominates over the single helical structure (whilst the opposite holds for the homochiral structures), which explains the aforementioned higher viscosity of the racemic bis-urea solution. This rationale constitutes a new basis to tune the macroscopic properties of the increasing number of supramolecular polymers reported to exhibit competing chiral nanostructures.

Introduction

Supramolecular polymers, i.e. chain-like assemblies of selfcomplementary monomers, are a class of compounds with unique properties and applications. [1][2][3][4] The directional but reversible noncovalent interactions responsible for their self-assembly provide a responsiveness [5] towards a range of stimuli that has contributed to the development of innovative catalysts, [6] optoelectronic materials, [7,8] self-healing materials, [9][10][11] chiroptical switches [12] and gels [13][14][15] amongst other applications. It was recently realized that the properties of these supramolecular polymers can be considerably diversified when their formation occurs in competition with other assemblies. Indeed, pathway complexity can lead to supramolecular polymorphism, when two supramolecular polymers form under kinetic [16][17][18][19][20][21] or thermodynamic [22][23][24][25][26][27][28][29][30][31][32][33][34][35] control. In the latter case, the transition between the two competing structures can be responsible for unusual behaviors such as self-assembly triggered by increasing the temperature, [28] or by dilution. [23] The modification of the supramolecular structure is particularly useful when it leads to a controllable and sudden change of macroscopic (such as rheological [22,30] or emissive [31,32] ) properties. In addition, the sharpness of the transition can actually be exploited to quantify weak intermolecular interactions [36,37] or solvation effects. [38] In this context, chirality is often used as a probe to allow monitoring subtle structural changes by circular dichroism (CD) spectroscopy. [16,[39][40][41][42][43] Moreover, chirality also directly impacts the thermodynamic stability of the assemblies. In the case of fibrillar self-assemblies or gels, heterochiral contacts can either increase or decrease the stability of the assemblies, compared to the homochiral situation. [44][45][46][47][48][49] In the case of supramolecular polymers made by stacking of a single row of monomers, rare cases exist where heterochiral dimers are more stable than homodimers, [50,51] but the usual situation is that heterochiral contacts lead to stereochemical defects. [52][53][54][55] When a heterochiral copolymer is formed by mixing two enantiomeric monomers, mismatches and helix reversals occur that destabilize the copolymer. If interpolymer interactions are negligible, the stability of the copolymer is therefore expected to increase with the enantiomeric purity as a result of the presence of fewer chiral defects. [56,57] Indeed, Shikata et al. have shown that an increase of the enantiomeric purity of a benzene-1,3,5-tricarboxamide monomer results in fewer defects along the supramolecular chains and thus in a higher viscosity of the solution. [58] Now, the situation is more complex when a chiral monomer can self-assemble through two distinct pathways: both resulting supramolecular polymers can potentially be affected by packing defects, so that the net result is difficult to predict. For example, we show herein that a racemic mixture of bis-urea stereoisomers forms thermodynamically stable heterochiral assemblies that result in a more viscous solution than that of homochiral assemblies formed by its chiral non-racemic stereoisomers. We show that this counter-intuitive result can actually be traced to a huge relative stability difference of the competing assemblies (ca 70°C), that can in turn be rationalized by a thermodynamic model that quantitatively accounts for the stereochemical defects.

Results and Discussion

The molecules described in this work contain four stereocenters (Fig. 1). In the following, we focus our attention only on the effect of the stereochemistry of the two valine units. The two ethylhexyl groups, that are used exclusively as racemates to improve solubility, [59] are assumed to have a negligible influence on the self-assembling bis-urea core. rac 2 -Val was synthesized by reacting an amine prepared from racemic valine with 2,4toluenediisocyanate and is therefore a mixture of four isomers (neglecting the isomerism due to the ethylhexyl groups). As previously reported, [60] it is readily soluble in non-polar solvents like toluene, where it forms a viscoelastic solution. Surprisingly, in the same conditions, SS-Val (derived from enantiopure S-valine) forms a fluid solution (Fig. 1). Small angle neutron scattering (SANS) measurements prove in both cases the formation of long and rigid rod-like assemblies (q -1 behavior at low wavevector q, Fig. 2). However, the scattered intensities differ by a factor of two at low q, indicating a difference in linear density of these one-dimensional supramolecular assemblies. A fit of the data with the form factor of long rigid rods with a circular cross section and a uniform scattering length density yields the values summarized in Table 1. These values indicate the presence of two (respectively one) monomers in the cross-section for rac 2 -Val (respectively SS-Val). The thicker rods can be expected to have a slower breaking dynamics, [22] which is qualitatively consistent with the obvious difference in viscosity of the samples (Fig. 1). SANS data show that at higher temperature (70°C), rac 2 -Val forms the same rods as SS-Val does at 20°C (Table 1). This actually suggests that RR, SS (and possibly RS and SR) monomers mix together and play a role in tuning the respective stability of the two rod-like structures. The influence of the optical purity of the monomers on the rheological properties of the solution is also observed in methylcyclohexane (MCH): mixing at room temperature the fluid solutions of SS-Val and RR-Val led to a viscous solution of rac-Val (see movie as supplementary material), which further demonstrates the dynamic nature of the supramolecular assemblies. This phenomenon was characterized more quantitatively by circular dichroism (CD) in MCH. Fig. 3a shows the CD spectrum of SS-Val in MCH at 20°C. The presence of a strong negative band at 233 nm is a clear sign of supramolecular chirality because this band is absent at high temperature (120°C). FTIR spectroscopy (Fig. S1) confirms that SS-Val forms hydrogen bonded assemblies in MCH at 20°C and SANS (Fig. S2, Table 1) proves that they consist in long rigid rods with a single monomer in the cross-section. Interestingly, a decrease in temperature triggers a significant change in the shape of the CD (Fig. 3a), UV-Vis (Fig. 3b) and FT-IR (Fig. S1a) spectra. These combined spectroscopic and scattering analyses prove that SS-Val selfassembles into helical rod-like objects, with two molecules (respectively, one molecule) in the cross-section at low (respectively, high) temperature. Their precise supramolecular structure has no consequence on the following analysis, but we surmise that SS-Val forms a double helix at low temperature and a single helix at high temperature, that are similar to the ones modelled for a structurally-related bis-urea monomer. [37] Figure 2. SANS analyses of bis-ureas rac 2 -Val and SS-Val performed at 20°C and at a concentration of ca. 6 g/L (10 mM) in toluene-d8. The curves are fitted according to the form factor for infinitely long [61,[START_REF]behaviour at low q is probably due to repulsive interactions between rods and precludes determination of their length[END_REF] rigid rods with a circular cross section and a uniform scattering length density. These features are reminiscent of the double helix and single helix structures found for similar bis-ureas (recalled as insets). [37] Table 1. Geometrical radius (r), linear density (nL) and number of bis-urea molecules (n) in the cross-section of the cylindrical objects deduced from a fit [61] of the SANS data at a concentration of ca. 6 g/L (10 mM). [a] Data from [60] . The transition between the two supramolecular structures of SS-Val is revealed by an abrupt change in the intensities of CD (Fig. 4a) and FTIR (Fig. S1b) signals in the range from ca. -10 to 0°C. The transition temperature (T**) has thus been measured as the midpoint of the jump in CD signal intensity for various mixtures of SS-Val and RR-Val (Fig. 4a). Since the racemic mixture is CD silent, the transition was measured by calorimetry (nDSC) (Fig. S3). [63] Fig. 4b shows that the CD, FTIR and nDSC data are consistent and that the evolution of the transition temperature is a strong function of the enantiomeric excess (ee). The racemic mixture (rac-Val) has a transition temperature that is 70°C higher than SS-Val and the relationship between T** and the ee is highly ) q (Å -1 )

SS-Val rac 2 -Val N N R N N R O H H
q -1 rac2-Val SS-Val
nonlinear. It confirms that SS-Val and RR-Val co-assemble and form heterochiral helical structures. Mixing RR-Val with SS-Val leads to a shift of the transition to higher temperatures: it therefore means that mixing the isomers either strongly stabilizes the low temperature structure or strongly destabilizes the high temperature structure. The stability of one-dimensional assemblies of chiral monomers is well-known to be affected by the presence of stereochemical defects such as mismatches (i.e. the presence of an enantiomer in a helix with its nonpreferred helicity) and helix reversals (i.e. a change of helicity within the assembly). [54] In the following part we propose a simple model to test if these intriguing results can be quantitatively accounted for. In the case of a scalemic mixture containing an arbitrary enantiomeric excess (ee) of SS-Val and RR-Val, the relative stability of the two kinds of supramolecular polymers is ruled by the free energy difference ∆𝐺 of binding a single monomer in the two types of assemblies [Eq. ( 1)]:

∆𝐺(𝑒𝑒) = 𝐺 '()' (𝑒𝑒) -𝐺 +,-(𝑒𝑒) (1)
We decompose each of these free energies into two contributions: the part that is due to stereochemical defects and the part that is insensitive to them, i.e. that does not depend on ee [Eq. ( 2)]: 𝐺 '()' (𝑒𝑒) = 𝐺 '()' ./0/123 (𝑒𝑒) -𝐺 '()'

4

(2)

the first contribution can be expressed as [Eq. ( 3)]: [START_REF] Van Gestel | [END_REF]65] 𝐺 '()' ./0/12 (𝑒𝑒) = ℎ '()' (𝑒𝑒) • 𝐻𝑅𝑃 '()' + 𝑚 '()' (𝑒𝑒) • 𝑀𝑀𝑃 '()' (3) where 𝐻𝑅𝑃 '()' (resp. 𝑀𝑀𝑃 '()' ) is the penalty for a helix reversal (resp. for a mismatch) in the high temperature supramolecular polymer, and ℎ '()' and 𝑚 '()' are the associated fractions of defects. Similarly, [Eq. ( 4)]: 𝐺 +,-(𝑒𝑒) = 𝐺 +,- ./0/123 (𝑒𝑒) -𝐺 +,- 4

with [Eq. ( 5)]:

𝐺 +,- ./0/123 (𝑒𝑒) = ℎ +,-(𝑒𝑒) • 𝐻𝑅𝑃 +,-+ 𝑚 +,-(𝑒𝑒) • 𝑀𝑀𝑃 +,- (5)
Therefore, [Eq. ( 6)]: ∆𝐺(𝑒𝑒) = ∆𝐺 ./0/123 (𝑒𝑒) -∆𝐺 4 (6) with [Eq. ( 7)]: ∆𝐺 ./0/123 (𝑒𝑒) = ℎ '()' (𝑒𝑒) • 𝐻𝑅𝑃 '()' + 𝑚 '()' (𝑒𝑒) • 𝑀𝑀𝑃 '()' -

ℎ +,-(𝑒𝑒) • 𝐻𝑅𝑃 +,--𝑚 +,-(𝑒𝑒) • 𝑀𝑀𝑃 +,- (7)
Since the fraction of defects is directly fixed by the value of the associated penalty (through the Boltzmann distribution), we see that the influence of ee can be quantitatively described by the four energetic parameters: 𝐻𝑅𝑃 '()' , 𝑀𝑀𝑃 '()' , 𝐻𝑅𝑃 +,-and 𝑀𝑀𝑃 +,-.

We then consider the free energy difference between the scalemic and the enantiopure systems [Eq. ( 8)]: ∆∆𝐺(𝑒𝑒) = ∆𝐺(𝑒𝑒) -∆𝐺(0) = >ℎ '()' (𝑒𝑒) -ℎ '()' (0)? • 𝐻𝑅𝑃 '()' + >𝑚 '()' (𝑒𝑒) -𝑚 '()' ( 0

)? • 𝑀𝑀𝑃 '()' -[ℎ +,-(𝑒𝑒) -ℎ +,-(0)] • 𝐻𝑅𝑃 +,--[𝑚 +,-(𝑒𝑒) -𝑚 +,-(0)] • 𝑀𝑀𝑃 +,- (8)
This free energy difference is directly related to the experimentally determined transition temperatures and can be calculated according to the supramolecular balance principle [Eq. ( 9)]: [36] [37] ∆∆𝐺(𝑒𝑒) = ∆𝐻 31B+/C(1 D𝑇 * * 31B+/C(1 -𝑇 * * GG H 𝑇 * * GG I

where both transition temperatures and the enthalpy have been measured (Fig. 4b andS3). Therefore, the energetic penalties are constrained by equation ( 10):

>ℎ '()' (𝑒𝑒) -ℎ '()' (0)? • 𝐻𝑅𝑃 '()' + >𝑚 '()' (𝑒𝑒) -𝑚 '()' (0)? • 𝑀𝑀𝑃 '()' -[ℎ +,-(𝑒𝑒) -ℎ +,-(0)] • 𝐻𝑅𝑃 +,--[𝑚 +,-(𝑒𝑒) - 𝑚 +,-(0)] • 𝑀𝑀𝑃 +,-= ∆𝐻 31B+/C(1 D𝑇 * * 31B+/C(1 -𝑇 * * GG H 𝑇 * * GG I ( 10
)
Moreover, following the van Gestel model for majority-rules experiments, [START_REF] Van Gestel | [END_REF] the helix reversal and mismatch penalties can be determined from the variation of the net helicity as a function of ee. Therefore, we selected temperatures below and above the lowest and highest T** values, respectively, to ensure that a single supramolecular polymer is present and we measured the CD spectra of the scalemic supramolecular polymers, both for the single helix supramolecular polymer (at 80°C) and the double helix supramolecular polymer (at -10°C) (Fig. S4). The deduced net helicity at 233 nm is plotted in Fig. 5a (single helix) and Fig. 5b (double helix) as a function of the ee of the SS-Val / RR-Val mixture. These plots show that the extent of net helicity amplification is strong for double helices, for which the net helicity reaches 70% for mixtures only slightly biased from the racemic mixture (20% ee), while it is modest for single helices, for which net helicity varies nearly linearly versus ee. Fig. 5c shows the thermodynamic data of Fig. 4b that is replotted according to equation (9). Finally, we used the equations of the van Gestel model (applied both to the single and double helices) together with equation (10) to simultaneously fit the data. The best fit is represented on Figs. 5a to 5c and the values obtained for the energetic penalties are collected in Table 2. The extent of net helicity amplification is reflected in the values of the energetic parameters deduced from the fit: the double helix is characterized by both a higher HRP and a lower MMP than the single helix. Both features are known to be necessary for amplification. [START_REF] Van Gestel | [END_REF]66] In addition, the larger HRP value for the double helix is expected, since a helix reversal should be more unfavorable in a thicker, probably less flexible, supramolecular polymer.

Table 2. Energetic parameters for stereochemical defects in SS-Val / RR-Val mixtures, extracted from the fits (Fig. 5). With these energetic parameters in hand, it is now possible to address the initially unanswered question whether the higher transition temperature of the racemic mixture is due to the stabilization of the low temperature supramolecular polymer or to the destabilization of the high temperature supramolecular polymer. Fig. 6a and6c show the fraction of defects that are present at temperatures above or below the transition temperature and Fig. 6b and6d show the associated energetic penalties derived from equations ( 3) and (5). For each supramolecular polymer, mismatches are much more numerous than helix reversals (i.e. 𝑚 '()' > ℎ '()' and 𝑚 +,-> ℎ +,-), but their energetic weight are similar (Fig. 6b and6d) because of the much higher penalty for a helix reversal than for a mismatch (i.e. HRPhigh > MMPhigh and HRPlow > MMPlow). As expected, for both supramolecular polymers the racemic composition contains more stereochemical defects that the non-racemic one, however the high temperature supramolecular polymer is much more destabilized by this effect than the low temperature one, because the difference between 𝐺 '()' ./0/123 (Fig. 6b) and 𝐺 +,- ./0/123 (Fig. 6d) is more than twice larger for the heterochiral than for the homochiral supramolecular polymer). This implies that for heterochiral assemblies, the double helical structure predominates over the single helical one, while the opposite holds for homochiral assemblies. Subtle energetic differences governing the stability of the two competing helical structures thus translate into a major change of the macroscopic properties of solutions of racemic versus enantiopure monomers.

Conclusions

Stereochemical defects are expected to destabilize a supramolecular polymer formed by a racemic mixture of monomers, compared to the enantiopure supramolecular polymer. However, when the monomers can assemble into two competing supramolecular polymers, the situation can be more complex. SS-Val self-assembles into two distinct supramolecular polymers; and the supramolecular polymer that is stable at low temperature is more viscous than the one stable at high temperature. When RR-Val is introduced in the system, it creates stereochemical defects that destabilize both supramolecular polymers. However, it is not possible, a priori, to predict how the relative stability of the two supramolecular polymers should evolve. The general thermodynamic model that we propose allows to rationalize this behavior based on four energetic penalties (for chiral mismatches and for helix reversals associated to both polymers). In the present experimental system, the high temperature supramolecular polymer is much more affected by stereochemical defects than the low temperature one, so that the relative stability of the low temperature supramolecular polymer is increased: its range of stability is increased by ca. 70°C. Consequently, at room temperature, the high temperature (low viscosity) supramolecular polymer is dominant in the case of SS-Val, but the low temperature (high viscosity) supramolecular polymer wins over in the case of the racemic mixture of monomers. This rationale should help understand the influence of chirality on the macroscopic properties of the increasing number of competing supramolecular chiral nanostructures reported in the literature.

Figure 1 .

 1 Figure 1. Structure of bis-ureas and photograph of their solutions at 10 mM (ca 6 g/L) in toluene (see also video as supplementary material).

  nL (Å -1 ) n[d] rac 2 -Val[a]

  [b] Measured by nDSC. [c] Measured by CD. [d] Number of bis-urea molecules in the cross-section, assuming a repeat distance of 4.6 Å.

rac 2 -

 2 Val = SS-Val + RR-Val + RS-Val + SR-Val

Figure 3 .

 3 Figure 3. CD (a) and UV (b) spectra of SS-Val measured at various temperatures in MCH (1 mM).

Figure 4 .

 4 Figure 4. (a) CD intensity at 233 nm versus temperature for SS-Val / RR-Val mixtures of various ee (100% ee corresponds to SS-Val). Conditions: 1 mM in MCH. The dashed vertical lines materialize the temperatures selected for net helicity measurements (Fig. 5). (b) Transition temperature derived from CD (a) or from FTIR (Figure S1) or nDSC (Figure S3) measurements.

4 Figure 5 .

 45 Figure 5. Majority rules experiments for SS-Val / RR-Val mixtures. (a) Net helicity of the high temperature supramolecular polymer (measured at 80°C). (b) Net helicity of the low temperature supramolecular polymer (measured at -10°C). Conditions: 1 mM in MCH, l = 233 nm. The CD spectra are shown in Fig. S4. (c) Relative stability of two supramolecular polymers. Continuous lines represent the fits with the parameters shown in Table2.

Figure 6 .

 6 Figure 6. Fraction of defects (at 80°C (a) and at -10°C (c)) and their energetic contribution (at 80°C (b) and at -10°C (d)) (calculated according to equations (3) and (5)), as a function of the enantiomeric excess of SS-Val / RR-Val mixtures.

rac 2 -

 2 Subtle energetic differences (chiral mismatches and helix reversals) govern the stability of two competing helical structures and translate into a major change of the macroscopic properties of solutions of racemic versus enantiopure monomers. = SS-Val + RR-Val + RS-Val + SR-Val Val = SS-Val + RR-Val + RS-Val + SR-Val viscous! fluid

Acknowledgements

This work was supported by the French Agence Nationale de la Recherche (project ANR-12-BS08-0019 BalanceSupra). Raphaëlle Pique, Jacques Jestin (LLB, Saclay) and Nicolas Vanthuyne (iSm2, Marseille) are acknowledged for assistance with the synthesis, SANS experiments, and chiral HPLC analyses, respectively.


for a clarinet model. In particular the dynamic oscillation threshold is predicted analytically when the blowing pressure is linearly increased. However, the mouth pressure cannot grow indefinitely. During a note attack, after an increasing phase, the musician stabilizes the mouth pressure. In the present work, the analytical prediction of the dynamic oscillation threshold is extended to a situations in which the mouth pressure approaches a steady state pressure according to an exponential time profile. The predictions still show a good agreement with simulation of the simple clarinet-model. This situation is compared in terms of dynamic oscillation bifurcation.

Introduction

One of the main sucesses of clarinet models, even if extremely simplified, is that they can predict ranges of parameters such as blowing pressure and lip force where the instrument produces a sound. The oscillation threshold, i. e. the minimum blowing pressure at which there can be a sustained oscillation has been extensively studied [START_REF] Kergomard | Calculation of the spectrum of self-sustained oscillators using a variable troncation method[END_REF][START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF]. The oscillation threshold can be measured by applying a constant blowing pressure, allowing enough time to let the system reach a permanent regime (either non-oscillating or strictly periodic), and repeating the procedure for other constant blowing pressures. This is a static or stationary view of the threshold.

A model based on an iterated map can be used to predict the asymptotic (or static) behavior (in particular the static oscillation threshold) of a simplified clarinet model as a function of a constant mouth pressure. This procedure avoids the phenomenon of bifurcation delay [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF], a shift of the pressure at which the oscillation starts when the pressure is gradually increased over time. In such a situation, the value of blowing pressure at which sound is observed is called dynamic oscillation threshold, in opposition to the static oscillation threshold. The phenomenon has been observed by Bergeot et al. in numerical simulations [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] and experiments [START_REF] Bergeot | Response of an artificially blown clarinet to different blowing pressure profiles[END_REF]. Using a simplified clarinet model the dynamic oscillation threshold has also been predicted [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF][START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF] when the mouth pressure is linearly increased.

In most realistic situations (for example during a note attack) the pressure stabilises at a target value as the oscillations grow to an audible level. The present paper shows how analytical results of [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF][START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF] can be extended to predict the dynamic oscillation threshold of the system for an archetypal mouth pressure shape that smoothly approaches the target value exponentially. The mathematical procedure and a comparison of theoretical results with numerical simulations are presented in section 3. Then, in section 4, a preliminary investigation on the influence of the mouth pressure shape on the onset of oscillations is performed. The clarinet model and major results from [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF][START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF] are first briefly recalled in section 2.

State of the art 2.1 Clarinet Model

This model divides the instrument into two elements: the exciter and the resonator. The exciter is represented by a nonlinear function F also called nonlinear characteristic of the exciter, relating the pressure applied to the reed p(t) to the flow u(t) through its opening. The resonator (the bore of the instrument) is described by its reflection function r(t). p and u are two non-dimensional state variables that are sufficient to describe the state of the instrument.

The solutions p(t) and u(t) depend on the control parameters: γ proportional to the mouth pressure P m according to

γ = P m P M = P m kH (1)
where P M = kH represents the pressure needed to close the reed entrance (also used to normalize the pressure p(t)), where 1/k is the acoustic compliance of the reed and H its distance to the lay at rest. The other parameter is ζ which is related to the opening of the embouchure at rest according to the formula

ζ = Z c U A /P M = Z c wH 2 ρP M . (2)
Here, Z c is the characteristic impedance at the input of the bore, w the effective width of the reed, and U A the maximum flow admitted by the reed valve. Biting harder the embouchure reduces the value of ζ. For most of the analysis below, this parameter is maintained at 0.5, but the analysis can easily be reproduced for other values of ζ.

The nonlinear characteristic is provided by the Bernoulli equation describing the flow in the reed channel plus some additional hypothesis on the turbulent mixing within the mouthpiece [START_REF] Hirschberg | A quasi-stationary model of air flow in the reed channel of single-reed woodwind instruments[END_REF][START_REF] Hirschberg | Aero-acoustics of wind instruments[END_REF]. The model is extremely simplified by considering a straight resonator in which the eventual losses are independent of frequency. In the current work, losses are ignored in all calculations. The reed is considered as an ideal spring [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF][START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF][START_REF] Ollivier | Idealized models of reed woodwinds. part 2 : On the stability of two-step oscillations[END_REF][START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF][START_REF] Chaigne | Instruments à anche[END_REF]. With these assumptions, the reflection function becomes a simple delay with sign inversion. Using the variables p + and p -(outgoing and incoming pressure waves respectively) instead of the variables p and u, the system can be simply described by an iterated map [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF]:

p + n = G p + n-1 , γ . (3)
Function G can be determined from the nonlinear characteristic F, which is done by Taillard [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF] for ζ < 1. This function depends on the control parameters γ and ζ. The time step n corresponds to the round trip time τ = 2l/c of the wave with velocity c along the resonator of length l.

Using the universal properties of iterated maps [START_REF] Feigenbaum | The universal metric properties of nonlinear transformations[END_REF], useful information about the instrument behavior can be drawn from the study of the iteration function. Most of these studies are done in the context of static bifurcation theory, which assumes that the control parameter γ is constant. For instance, it is possible to determine the steady state of the system as a function of the parameter γ, and to plot a bifurcation diagram shown in red in Fig. 1, in terms of variable p + . When no losses are considered, we have γ st = 1/3. For all values of the control parameter γ below γ st the series p + n converges to a single value p + * corresponding to the fixed point of the function G, i.e. the solution of p + * = G (p + * ). When the control parameter γ exceeds γ st the fixed point of G becomes unstable and the steady state becomes a 2-valued oscillating regime.

Slowly linear time-varying mouth pressure 2.2.1 Dynamic bifurcation

A control parameter γ increasing linearly with time is taken into account by replacing Eq. ( 3) by Eqs. (4a) and (4b):

     p + n = G p + n-1 , γ n (4a) γ n = γ n-1 + . ( 4b
)
γ is assumed to increase slowly, hence is considered arbitrarily small [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]. When the series p + n is plotted with respect to parameter γ n the resulting curve can be interpreted as a dynamic bifurcation diagram and it can be compared to the static bifurcation diagram (Fig. 1).

Because of the time variation of γ, the system in Eqs. ( 4) is subject to the phenomenon of bifurcation delay [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Fruchard | Sur le retard à la bifurcation[END_REF]: the bifurcation point (in this case the oscillation threshold) is shifted from the static oscillation threshold γ st [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] to the dynamic oscillation threshold γ dt [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]. The difference γ dt -γ st is called the bifurcation delay. In a previous work [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF], the bifurcation delay was found to depend very strongly on noise, in particular due to round-off errors made by the computer in numerical simulations, even if high precisions were used.

According to dynamic bifurcation theory, two operative regimes must be distinguished [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF]:

• The Deterministic Regime (DReg.): In this case, the noise does not affect the bifurcation delay which does not depend on the slope of the blowing pressure.

• The Sweep-Dominant Regime (SDReg.): In this case, the bifurcation delay is affected by the noise, becoming larger as the blowing pressure γ is increased quicker.

Articles [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF][START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF] provide an analytical study of the dynamic bifurcation of the clarinet model (i.e. Syst. ( 4)) based on a generic method given by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF]. The main results of these studies are theoretical estimations of the dynamic oscillation threshold of the clarinet: one for the DReg. [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] and one for the SDReg. [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF]. These expression are recalled below.

Dynamic oscillation threshold for the deterministic regime

A possible theoretical estimation of the dynamic oscillation threshold consists in identifying the value of γ for which the orbit of the series p + n escapes from a neighborhood The phenomenon of bifurcation delay is highlighted.

of arbitrary distance of an invariant curve φ(γ, ). More precisely, the dynamic oscillation threshold is reached when the distance between the orbit and the invariant curve becomes equal to . The invariant curve (i.e. invariant under the mapping (4), described for example in [START_REF] Fruchard | Bifurcation delay and difference equations[END_REF]) can be seen as the equivalent of a fixed point in static regimes. It satisfies the following equation:

φ (γ) = G (φ (γ -), γ) . (5)
Although Eq. ( 5) usually leads to mathematical expressions that cannot be calculated analytically, the invariant curve can be determined approximately with a perturbation method given by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF], leading to the following general form:

φ (γ) = n i=0 i φ i (γ) + o( n+1 ), (6)
where the zeroth order term of the series is the fixed point curve of the function G, φ 0 (γ) = p + * (γ). The procedure to obtain the theoretical estimation γ th dt of the dynamic oscillation threshold is as follows: a theoretical expression of the invariant curve is found for a particular (small) value of the increase rate (i.e. 1). Equations ( 4) are then expanded into a first-order Taylor series around the invariant curve and the resulting linear system is solved analytically. Finally, γ th dt is derived from the analytic expression of the orbit.

The analytic estimation of the dynamic oscillation threshold γ th dt is defined in [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]:

γ th dt + γ 0 + ln ∂ x G φ (γ -), γ dγ = 0, ( 7
)
where γ 0 is the initial value of γ (i.e. the starting value of the linear ramp).

Dynamic oscillation threshold for the sweepdominant regime

The effect of the noise can be taken into account by introducing an uniformly distributed random variable in the system described by Eqs. ( 4). This random variable is an additive white noise with an expected value of zero and variance σ 2 . In the case of a numerical simulation performed with finite precision,1 σ = 10 -precision .

The method to obtain the theoretical estimation of the dynamic oscillation threshold for the SDReg. (noted γth dt ) is detailed in [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF].

Because of noise, the bifurcation delay is reduced. Therefore, the main approximation of the method is to assume that the dynamic oscillation threshold is close to the static oscillation threshold γ st . Using this approximation, the expression of γth dt is:

γth dt = γ st + - 2 K ln π K 1/4 σ 5/4 , (8)
which is the theoretical estimation of the dynamic oscillation threshold for the SDReg., K is a constant that depends on the slope of ∂ x G(p + (γ), γ), the derivative of the iteration function at the fixed point.

Exponential variation of the mouth pressure

During a note attack, the mouth pressure cannot grow indefinitely, being stabilized to a targeted value before the oscillations grow to an audible level. This section is devoted to show how results presented in sections 2.2.2 and 2.2.3 can be extended to predict the dynamic oscillation threshold of the system for a profile in which the parameter approaches asymptotically a target value through an exponential function.

Prediction of the dynamic oscillation threshold: mathematical procedure

When the mouth pressure follows an exponential function, the system is described by the following system of difference equations:

     p + n = G p + n-1 , γ n (9a) γ n = a γ n-1 + γ M (1 -a), (9b)
where γ M is the targeted mouth pressure (it is always equal to 1 in this work). Eq. (9b) describes the exponential variation of the mouth pressure, with γ 0 = 0 and noting =ln[a], its formal solution is given by:

γ n = γ M 1 -e -n , (10)
Fig. 2 shows an example of numerical simulation performed on the system described by Eqs. [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF].

Using a change of parameters, Eqs. ( 9) can take the form of (4). The iterative function G is replaced by H:

H (x, η) = G (x, γ(η)) , (11)
with a linearly increasing parameter η: The new system is written:

η (γ) = ln γ M γ M -γ =⇒ γ (η) = γ M 1 -e -η . (12
     p + n = H p + n-1 , η n (13a) η n = η n-1 + . ( 13b
)
Using Eqs. ( 7) and ( 8), it possible to predict the dynamic oscillation thresholds, noted η th dt (for DReg.) and ηth dt (for SDReg.), of the system of Eqs. [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF]. These thresholds are then expressed in terms of mouth pressure using Eq. ( 12

Benchmark of theoretical estimators for the dynamic threshold

In this section, the above theoretical predictions of the dynamic threshold obtained for an exponential variation of the mouth pressure are compared to numerical simulations. A numerical dynamic threshold η num dt is estimated as the The method presented in this article provides an extension of the estimation of dynamic thresholds for exponentially increasing parameters. The present method can, in principle, be used for any other profile of timevarying parameter that can be described analytically.

Previous works have shown that, for sufficient precision, the dynamic threshold is independent of the rate of variation of the parameter (Deterministic regimes). A quick extrapolation of these results might have led to the conclusion that in a generic profile, even if there is a variation of the rate of change of γ, there would be no significant changes in the dynamic threshold. However the results in this article show that this is not the case, as the deterministic regime has a threshold (approximately 0.7) that is smaller than that of the linearly increasing profile (approximately 0.9).

In real cases, however, the system is always far from a deterministic regime, as the numerical or turbulence noise introduces a stochastic variation in the parameter then brings the system into the sweep-dominant regime. In these cases (similarly to the deterministic regime), exponential shape appears to provide faster onset of oscillations. Obviously, additional works must be performed to state definitive conclusions about the influence of the mouth pressure shape on the onset of oscillations in a clarinet.

Figure 1 :

 1 Figure 1: Comparison between static and dynamic bifurcation diagram as functions of γ n . = 2 • 10 -3 , ζ = 0.5.The phenomenon of bifurcation delay is highlighted.

Figure 2 :

 2 Figure 2: Numerical simulation performed on the Syst. (9).The parameters used are: γ M = 1, a = 0.995 ( = 0.005), ζ = 0.5, γ 0 = 0 and p + 0 = G (0; γ 0 ).

Table 1 :

 1 Notation for thresholds of oscillation.

	):

The precision in here referred as the number of decimal digits used by the computer.

value for which the distance between the simulated orbit of Syst. [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF] and its invariant curve is first larger than . The comparison is carried out as a function of parameter . Results are shown in terms of mouth pressure γ in Fig. 3 with Γ num dt = γ η num dt and for several values of the precision. The choice of the precision is possible using mpmath, an arbitrary precision library for Python.

Similarly to the case of a linear variation of the mouth pressure γ [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF], the evolution of Γ num dt with respect to allows to distinguish the deterministic and sweep-dominant regimes: in certain areas of the figures, the dynamic bifurcation threshold does not depend strongly on , this is the DReg., while in other areas the dynamic bifurcation threshold depends on , this is the SDReg.. The lower is the precision, the larger is the value of for which the DReg. appears. For each regime the theoretical results Γ th dt and Γth dt provides a good estimation of the dynamic oscillation threshold of the clarinet model ( 9).

Linear vs. exponential variation of the mouth pressure

In this section, numerical dynamic threshold obtained for an exponential variation of the mouth pressure (Γ num dt ) is compared to numerical dynamic threshold obtained for a linear variation of the mouth pressure (γ num dt ). One can raise the question of the chosen parameter to carry out this comparison. As a preliminary investigation, we choose here the time (noted N) needed to reach 99 per cent of the target value γ M (see Fig 4).

For a linear variation of the mouth pressure with an increase rate , we have:

Inverting Eq. ( 10) gives:

Therefore, for an exponential variation of the mouth pressure, N is defined by:

To illustrate the definition of the parameter N, Fig. 4 shows a linear and an exponential functions plotted with the same value of N. The comparison is depicted in Fig. 5. Fig. 5(a) compares the dynamic oscillation thresholds Γ num dt and γ num dt : for the DReg., γ num dt is larger than Γ num dt while, for SDReg., the opposite is noticed.

To complete the study, it is also interesting to compare the times N num dt and n num dt (computed through Eq. ( 15)) needed to reach Γ num dt and γ num dt (see Fig. 5(b)). For DReg. as well as for SDReg., exponential shape appears to provide faster onset of oscillations.