Julien Pouysegur
email: julien.pouysegur@institutoptique.fr

Florent Guichard

Yoann Zaouter

Marc Hanna

Frédéric Druon

Clemens Hönninger

Eric Mottay

Patrick Georges

Hybrid high-energy high-power pulsewidthtunable picosecond source

A hybrid ytterbium-doped fiber-bulk laser source allowing the generation of 3 ps, 350 μJ, 116 MW peak power Fourier transform-limited pulses at 50 kHz repetition rate and 1030 nm wavelength is described. Pulse duration tunability is provided by an adjustable spectral compression-based seeder system. Energy scaling capabilities of the architecture by use of the divided-pulse amplification method are investigated. This source provides a robust, compact, and versatile solution for applications such as RF photocathode guns, x-and γ-ray generation by inverse Compton scattering, and optical parametric chirped-pulse amplification pumping.

High-repetition rate, high-energy picosecond sources are required for the development of secondary sources, such as photocathode RF guns, x-and γ-ray generation using inverse Compton scattering, and optical parametric chirped-pulse amplification (OPCPA) systems. These applications have noteworthy operating points in the 1-20 ps laser pulse duration regime, and often specify a narrow spectral bandwidth. For electron gun and OPCPA applications [1][2][3], the narrow spectral width is dictated by the need to use the second, third, or fourth harmonic of the laser initial frequency. Nonlinear frequency conversion is made considerably simpler if the initial pulses are close to Fourier transform-limited (FTL), relaxing constraints on the spectral acceptance of the process. X-ray sources using inverse Compton scattering use resonant cavity enhancement of the laser radiation, for which a narrower spectrum allows much larger field enhancement factors [4,5].

A widely adopted solution to amplify picosecond pulses relies on the use of Nd-doped materials, such as Nd:YVO 4 . Such amplifiers benefit from a large achievable gain and moderate setup complexity [6][7][8]. By cascading several amplification stages or using multipass geometries, energies of several hundreds of microjoules can be obtained in 10 ps pulses [9,10]. However, the narrow spectral emission cross section of neodymium does not allow the generation of pulse durations below 10 ps [10,11]. Moreover, average power scaling of these systems is limited by the relatively large quantum defect.

In contrast, systems based on ytterbium-doped materials allow average power scaling due to the low quantum defect of the pumping scheme, and broader spectral emission cross sections permit the generation of much shorter pulses, well in the femtosecond-picosecond regime [12,13]. For example, a system delivering 1.3 mJ, 7 ps pulses at 800 kHz repetition rate, corresponding to an average power beyond 1 kW, has been demonstrated [13], although in a relatively complex setup.

Ytterbium-doped fibers amplifiers (YDFAs) can also be used to amplify picosecond pulses, with the associated benefits of large gain, beam quality, and average power scaling, but the small core dimensions restrict the pulse energy below 100 μJ [14,[START_REF] Zaouter | Proc. SPIE[END_REF]. The chirped-pulse amplification technique (CPA) that is routinely used to scale the energy of femtosecond systems cannot be easily implemented for FTL picosecond pulses due to the narrow spectral bandwidth. However, an equivalent attractive alternative, denoted as divided-pulse amplification (DPA), has recently been used in this context [START_REF] Zhou | [END_REF]17], and in high-energy femtosecond fiber amplifiers [18].

In this Letter, we demonstrate a FTL picosecond source composed of a seeder system tunable in pulse duration and a hybrid fiber-bulk amplifier able to provide large gain, high pulse energy, and high average power. Pulse-width tunability is obtained by exploiting the spectral compression effect induced by self-phase modulation (SPM) in an optical fiber. The DPA architecture is used to scale the pulse energy beyond the self-focusing threshold both in the fiber amplifier and in the Yb:YAG bulk amplifier. This system has the advantage to be relatively simple compared to already existing sources used for the same applications [11,19]. This results in a source capable of delivering 3 ps, 350 μJ pulses at a repetition rate of 50 kHz, corresponding to 17.5 W average power. These performances are obtained in a compact and robust setup, and can be easily adjusted in the 3-20 ps range of pulse duration, in energy (10-350 μJ), and in average power (10-50 W), making it a particularly versatile source.

The experimental setup, schematically shown in Fig. 1, consists of a front-end generating tunable-in-duration picosecond pulses followed by a fiber power amplifier and an Yb:YAG booster for the generation of energetic output pulses. The front-end is composed of an Yb:KYW bulk oscillator (Mikan, Amplitude Systemes) operating at 54 MHz repetition rate, and delivering 250 fs pulses centered at 1030 nm with 1.5 W average power. Pulses are then temporally stretched to picosecond durations using a 1750 l/mm grating-based arrangement introducing negative chirp. A fiber-based acoustooptic modulator is used to decrease the repetition rate of the oscillator before amplification in a core-pumped YDFA using a 6 μm core and a 2 m long fiber. The repetition rate of the fiber amplifier is 4 MHz. When negatively chirped pulses are subjected to SPM in this amplifier, spectral components at both ends of the spectrum are redistributed toward the central wavelength, leading to a spectral narrowing of the pulses [20]. After a propagation length that depends on the initial chirp and strength of the nonlinearity, FTL pulses can be obtained. By controlling the input chirp of the pulse and gain of the amplifier, FTL pulses with continuously tunable pulse width in the range 3-30 ps are easily generated. Figure 2 shows the oscillator and picosecond pulse spectrum after spectral compression-for initial group delay dispersion of -0.48 ps 2 and -2.25 ps 2corresponding to a pulse duration of 3 ps and 20 ps, respectively. The obtained full width at half-maximum spectral bandwidths are 0.7 and 0.19 nm, respectively, confirming the quasi-FTL nature of these pulses. The typical pulse energy out of this first pulse-shaping amplifier stage is 5-30 nJ.

Since the generation of FTL picosecond pulses relies on a nonlinear process, a legitimate concern arises over the stability of the narrow spectrum. To investigate this aspect, a long-term measurement of the output spectrum is carried out in the configuration where the amount of SPM is the largest, corresponding to long pulse durations. Figure 3 displays the spectrum of 20 ps pulses and the average power at the output of the first YDFA stage recorded over 9 h, showing power fluctuations of 0.3% RMS and spectral width fluctuations of 3% RMS. The use of a highly stable bulk oscillator together with wellintegrated fiber components allows this stability, and residual fluctuations are attributed to beam pointing fluctuations at the input of the fiber subsystem.

Further described experiments are carried out with 3 ps pulses. The first power amplifier is based on a 0.8 m long 85/265 μm double-clad distributed-mode filtering rod-type ytterbium-doped fiber pumped with a 50 W fiber-coupled high-power pump diode system at 976 nm.

It is operated in a double-pass configuration. Despite the very large mode area of the fiber, SPM can induce spectral broadening in this stage depending on the output energy. A convenient way to mitigate this effect consists of partially distributing the spectral compression process in the power amplifier by adjusting the gain of the shaping stage. This allows a maximum pulse energy of 1.5 μJ for 3 ps pulses, corresponding to 500 kW peak power, while preserving FTL pulses at the output.

Further energy scaling of this stage is investigated based on the use of theDPA concept. This technique relies on the division of the input pulse into several less energetic temporal replicas [13]. After amplification, these replicas are combined to regenerate a single output pulse. The replicas for this power amplifier stage are generated with a set of three properly oriented yttrium vanadate crystals, which exhibit a large group delay birefringence of 0.7 ps∕mm. The length of the smallest crystal is 10 mm corresponding to a differential group delay of 7 ps, twice the input pulse duration. The two other crystal lengths are 20 and 40 mm (by combining two 20 mm crystals), leading to the generation of eight pulse replicas. In each crystal, the incident polarization state is linear and oriented at 45°o f the optical axes, generating at the output two temporally separated sets of replicas of equal amplitude and orthogonal polarization. Passive recombination is achieved through a 90°polarization rotation (using a Faraday rotator), which ensures a complementary propagation path inside the birefringent crystals on the way back, and thus a complete compensation of the previously induced delays. Nonperfect recombination can originate in differential effects among replicas such as a phase shift or gain. In this case, part of the energy of the replicas is directed toward the orthogonal output polarization state, denoted as the uncombined output. Parasitic pulses in the time domain can also appear on the combined output. Energy scaling results using DPA are shown in Fig. 4. At low energy per replica and for two replicas, the combining efficiency, defined as the combined power divided by total output power, is 99%. This efficiency decreases to 96% for eight temporal replicas, and we attribute this decrease to slight misalignment of crystal orientations in the dividing/combining stages. At the maximum energy of 1.5 μJ per replica, the efficiency decreases from 96% for two replicas down to 88% for eight replicas. This is attributed to the increased sensitivity of the system to nonperfect replica amplitude equalization that can be turned into phase discrepancies among replicas through accumulated SPM in the amplifier, as detailed in [21]. However, combined output energy of 9.3 μJ is achieved for 3 ps pulses, corresponding to 3 MW peak power, a record value for a quasi-FTL ps fiber amplifier [14,[START_REF] Zaouter | Proc. SPIE[END_REF]17] To further scale the energy, an additional amplifier based on a 10 mm long Yb:YAG crystal is used. Because of the large gain provided by this bulk booster stage, the DPA divider/combiner is used at the input of the bulk booster. Indeed, the output energy provided by the power fiber amplifier with DPA-using one vanadate crystal-generating two replicas, approximately 3 μJ, is sufficient to ensure efficient energy extraction in the bulk booster, and DPA is needed at this stage to provide energy scaling capabilities. The DPA divider for the bulk Yb:YAG booster generates four replicas to accommodate higher peak power: the first temporal division is based on a combination of thin-film polarizer and free-space delay line (10 cm∕330 ps), and the second division is based on a 20 mm long yttrium vanadate crystal (15 ps). This allows the final combining process to take place in free space, thereby avoiding self-focusing issues.

The input beam is focused inside the Yb:YAG crystal down to a 400 μm waist diameter. The crystal is pumped by a diode system emitting 170 W of average power at 969 nm, and coupled to a 200 μm, 0.22 NA pump fiber. The output facet of the pump fiber is imaged inside the crystal with a diameter of 500 μm. Heat removal in the crystal is operated by water cooling of the top and bottom facets of the crystal through a copper mount. The crystal is thermally contacted to the copper mount by a sheet of indium.

We first investigate the output power/energy of the bulk amplifier as a function of repetition rate at a maximum pump power of 170 W, as shown in Fig. 5. An additional acoustooptic modulator located between the rod-type amplifier and the bulk booster allows us to make this systematic study. The output pulse energy increases from 45 μJ at 1 MHz repetition rate to 350 μJ at 50 kHz, corresponding to average powers of 45 and 18 W, respectively. The 350 μJ energy corresponds to a peak power of 116 MW, and the Kerr focal length associated with this peak power is <1 cm in YAG. It is clear that severe self-focusing would occur without using the DPA system. With the DPA activated, the B-integral is estimated around 0.22 rad∕replica. In order to decouple, identify, and compensate the thermal lens and the Kerr lens, the setup is aligned in two steps. First, we operate at full pump power but low seed energy to observe the thermal lens only. The backreflecting mirror position is then optimized to compensate it. Afterward, nonlinear effects can be observed by increasing the input signal. The beam propagation changes imparted by the Kerr lens during first and second passes and is more difficult to fully compensate.

The combined power, uncombined power, and combining efficiency are plotted as a function of pump power at 50 kHz repetition rate in Fig. 6. The DPA combining efficiency of the Yb:YAG booster is above 97% in all cases, because of the low value of B-integral estimated to be 0.3 rad per replica (350 μJ output pulse energy). It is also worth noticing that all depolarization effects are also low since they are included in the 3% combination losses. This illustrates the relevance of DPA for bulk amplifiers pumped with high-brightness diodes, allowing energy-scaling beyond single-pulse self-focusing and optical damage limitations. No noticeable evolution is observed between input and output spectra of the booster amplifier, as shown in Fig. 5. Indeed, spectral evolution essentially occurs during spectral compression, which takes place in the fiber system. The amplified spontaneous emission (ASE) appears to be negligible, as shown in Fig. 5 representing the amplified spectrum in log-scale.

Finally, we investigate spatial properties of the output pulses at 50 kHz repetition rate and maximum energy of 350 μJ. The spatial quality is also pristine, with M 2 values below 1.15 in both directions and a Gaussian beam profile. For a higher repetition rate, above 40 W output power, spatial distortions start to appear due to additional heat load induced by strong pump absorption and high gain saturation. The thermal lens limitation is not related to its dioptric power value but to aberrations that cause beam profile degradation.

We also investigated the energy extraction with 20 ps pulses (see Fig. 7). Thanks to the reduced peak intensity we were able to extract 720 μJ from the Yb:YAG booster without implementing DPA. Damage threshold was reached slightly above 860 μJ for a single pulse corresponding to a fluence around 500 mJ∕cm 2 . The use of DPA allows us to operate significantly below the damage threshold.

In conclusion, we report a hybrid high-power versatile source of FTL picosecond pulses based on a rod-type fiber amplifier followed by a bulk Yb:YAG booster. Energy scaling using the DPA technique for the fiber amplifier allows obtaining record performances for this type of system. This method is also advantageously used to scale the energy of the bulk booster, allowing the generation of 350 μJ, 3 ps pulses for a peak power of 116 MW. Higher energy has been reached up to 720 μJ with 20 ps pulse duration for a peak power of 36 MW without any divided pulse scheme. We believe that this type of versatile architecture will allow further energy scaling to the millijoule level mostly by adding division stages, since 32 replica systems have already been demonstrated [START_REF] Zhou | [END_REF]. It is worth noticing that with such a hybrid configuration based on a fiber-bulk amplifier, we measured in hybrid industrialized systems power fluctuation of 0.3% rms and beam pointing stability of 30 μrad RMS over 15 h. The bulk booster allows simple energy scaling of over an order of magnitude at the output of a CPA free picosecond fiber system [11] while preserving a very clean and narrow spectrum compatible with highly efficiency frequency conversion.

Funding. French National Research Agency (ANR) (Femtocryble program).

Fig. 1 .

 1 Fig. 1. Schematic laser setup. HWP, half-wave plate; FR, Faraday rotator.

Fig. 2 .

 2 Fig. 2. (a) Initial oscillator spectrum (blue), compressed spectrum for 3 ps (red) and 20 ps (green) pulses. (b) Together with corresponding autocorrelation traces.

Fig. 3 .

 3 Fig. 3. Long-term stability of the front-end for 20 ps pulse generation. Spectral shape at sampled times (right); average power stability of the oscillator (top).

Fig. 4 .

 4 Fig. 4. Evolution of (a) the combined output energy and (b) combining efficiency as a function of energy per replica for different numbers of replicas in the fiber power amplifier.

Fig. 5 .

 5 Fig. 5. (a) Output spectrum at 50 kHz, (b) output energy and output power as a function of repetition rate at the output of the bulk booster stage for 3 ps pulse duration. Inset: far-field beam profile at 350 μJ.

Fig. 6 .

 6 Fig. 6. Combined, uncombined output power, and combining efficiency as a function of pump power at 50 kHz at the output of the bulk booster.

Fig. 7 .

 7 Fig. 7. Output energy and output power as a function of repetition rate at the output of the bulk booster stage for 20 ps pulse duration. Inset: far-field beam profile at 720 μJ.


Landrieu Loic

Guillaume Obozinski

Cut Pursuit: fast algorithms to learn piecewise constant functions

We propose working-set/greedy algorithms to efficiently find the solutions to convex optimization problems penalized respectively by the total variation and the Mumford Shah boundary size. Our algorithms exploit the piecewise constant structure of the level-sets of the solutions by recursively splitting them using graph cuts. We obtain significant speed up on images that can be approximated with few level-sets compared to state-of-the-art algorithms.

Introduction

Estimation or approximation with piecewise constant functions has many applications in image and signal processing, machine learning and statistics. In particular, the assumption that natural images are well modeled by functions with bounded total variation motivates the use of the latter as a regularizer, which leads to piecewise constant images for discrete approximations. Moreover a number of models used in medical imaging (El-Zehiry and Elmaghraby, 2007) assume more directly piecewise constant images. More generally piecewise constant models can be used for compression, for their interpretability and finally because they are typically adaptive to the local regularity of the function approximated [START_REF] Wang | Trend filtering on graphs[END_REF]. Piecewise constant functions display a form of structured sparsity since their gradient is sparse.

Both convex and non-convex formulations have been proposed to learn function with sparse gradients, the most famous being (a) the formulation of [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], hereafter referred to as ROF, who proposed to regularize with the total variation, and (b) the formulation of [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] who proposed to penalize with the total length of discontinuities of piecewise smooth functions. A fairly large literature is devoted to these formulations mainly in the image processing and optimization literature. Although the connection between the total variation, the Mumford-Shah energy and graph cuts is today well established, algorithms that leverage this connection are relatively recent. In particular, for ROF, [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF]; [START_REF] Goldfarb | Parametric maximum flow algorithms for fast total variation minimization[END_REF] use the fact that the problem can be formulated as a parametric max-flow. El-Zehiry and Grady (2011) use graph cuts to solve the formulation of Mumford and Shah for the case of two components.

The literature on sparsity in computational statistics and machine learning has shown how the sparsity of the solutions sought can be exploited to design algorithms which use parsimonious computations to solve the corresponding large scale optimization problem with significant speed-ups [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF]. Our work is motivated by the fact that this has to the best of our knowledge not been fully leveraged to estimate and optimize with piecewise constant functions. In the convex case, the algorithm proposed to exploit sparsity are working set algorithms and the related Frank-Wolfe. In the non-convex case forward selection algorithms such as OMP, OLS, CoSamp, FoBa and others have been proposed [START_REF] Mallat | Adaptive timefrequency decomposition with matching pursuits[END_REF][START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF][START_REF] Zhang | Adaptive forward-backward greedy algorithm for sparse learning with linear models[END_REF] 1 .

It is well understood that algorithms for the convex and non-convex case are in fact fairly related. In particular, for a given type of sparsity the forward step of working set methods, Frank-Wolfe or greedy algorithms is typically the same, and followed by the resolution of a reduced problem.

Given their similarity, we explore in this paper both greedy and working set strategies to solve respectively the problems penalized by the Mumford-Shah penalty for piecewise constant functions (aka minimal partition problem) and the optimization problem regularized by the total variation. In the convex case, our algorithms do not apply only to the case where the data fitting term is the MSE or a separable smooth convex function, for which some efficient algorithms implicitly exploiting sparsity exist [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF][START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF][START_REF] Kumar | Active-set methods for submodular optimization[END_REF], but also to a general smooth convex term.

Our algorithms are very competitive for deblurring and are applicable for the estimation of piecewise constant functions on general weighted graphs.

Notations

We will denote G = (V, E, w) an unoriented weighted graph whose edge set is of cardinality

m and V = [1, • • • , n].
For convenience of notation and proofs we encode the undirected graph G as a directed graph with oriented edges in both directions in place of unoriented edges. For a set of nodes A ⊂ V we denote 1 A the vector of {0, 1} n such that [1 A ] i = 1 if and only if i ∈ A. For F ⊂ E a subset of edges we denote w(F ) = (i,j)∈F w ij .By extension for two subsets A and B of V we denote w(A, B) = w(A × B ∩ E) the weight of the boundary between those two subsets. Finally we denote C 2 V the set of all partition of V into connected components in the graph G.

General problem considered

Problem formulation

We consider in this work the problem of minimizing functions Q of the form Q(x) := f (x) + λΦ(x) with f differentiable and Φ(x) a penalty function that decomposes as : Φ(x) = ij∈E w ij φ(x i -x j ) with φ : R → R + a sparsity inducing function such that φ(0) = 0. The general problem writes

min x ∈ R n f (x) + λ 2 (i,j)∈E w ij φ(x i -x j ).
(1)

The function φ is typically the absolute value, which corresponds to the total variation, or one minus the Kronecker delta at 0, which leads to the Mumford-Shah penalty for piecewise constant functions. More generally, for functions φ that have a nondifferentiability at 0, the solution of (1) has sparse gradient and is thus constant on the elements of a certain coarse partition of V . We therefore reformulate the problem for candidate solution that have that property.

Decomposition on a partition

Any x ∈ R n can be written as

x = k i=1 c i 1 Ai with Π = {A 1 , • • • , A k } ∈ C a partition of V into k connected components and c ∈ R k . Con- versely we say that x can be expressed by parti- tion Π = (A 1 , • • • , A k ) if it is in the set span(Π) := span(1 A1 , • • • , 1 A k ) = { k i=1 c i 1 Ai | c ∈ R k }.
We denote x Π = arg min z∈span(Π) Q(z) the solution of (1) when x is constrained to be in span(Π). With this notation, we can rewrite problem (1) as the problem of finding an optimal partition Π :

Π = arg min Π∈C Q(x Π ) (2)
We say that a partition Π is coarse if its cardinality k is such that k n.

Before presenting our approach, we review some of the main ideas relevant in the related literature.

1.3 Related work [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] describe an image as simple if it can be expressed as a piecewise-smooth function, i.e. if the domain of the image can be partitioned into regions with short contours and such that the image varies smoothly inside of each region. The socalled minimal partition problem refers to the particular case where the functions are taken to be constant inside of each region. The setting in which the number of regions is known prior to the optimization, typically two, is known as the Chan-Vese problem and was first solved using active contour methods [START_REF] Kass | Snakes: Active contour models[END_REF]. [START_REF] Chan | Active contours without edges[END_REF] propose a level-set based method for the binary case, which has the advantage of foregoing edges and gradient completely, as they are typically very sensitive to noise. This method has since been extended to the so called multiphase setting where the number of phases, that is of level-sets of the function is a power of two [START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF]. The resolution of those problems was considerably sped up by the introduction of graph-cut methods, for binary phase [START_REF] El-Zehiry | Brain MRI tissue classification using graph cut optimization of the Mumford-Shah functional[END_REF] and in the multiphase setting (El-Zehiry and Grady, 2011).

A formulation which can be considered a posteriori very related to the one of Mumford and Shah was proposed by [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and consists in regularizing the estimated image with the total variation. Their celebrated approach boasts numerous applications in various fields such as vision and signal processing, and is solved in the recent literature using proximal splitting [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Raguet | A generalized forward-backward splitting[END_REF].

The relation between graph-cuts and the total variation goes back to [START_REF] Picard | Minimum cuts and related problems[END_REF] but has been only fully exploited recently, when [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF] and [START_REF] Goldfarb | Parametric maximum flow algorithms for fast total variation minimization[END_REF] among others, have exploited the fact that the ROF model can be reformulated as a parametric maximum flow problem; the latter entails that it can be solved by a divide-and-conquer strategy which requires to solve a sequence of max-flow problems on the same graph, thus allowing for efficient reuse of partial computation with a push-relabel algorithm. These results on the total variation are actually an instance of results that apply more generally to submodular functions [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF]. Indeed, the intimate relation existing between the total variation and graph-cuts is due fundamentally to the fact that the total variation is the Lovász extension of the value of the cut which is a submodular function. Recent progress made on efficient optimization of submodular function produced new fast algorithm for the total variation [START_REF] Kumar | Active-set methods for submodular optimization[END_REF][START_REF] Jegelka | Reflection methods for user-friendly submodular optimization[END_REF].

Problems regularized by the total variation or the Mumford-Shah energy are both related to Potts model. Indeed, if the values of the level-sets are quantized, the corresponding energy to minimize is that of a discrete valued conditional random field (CRF), with as many values as there are quantization levels [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF][START_REF] Chambolle | An introduction to total variation for image analysis[END_REF]. A number of optimization techniques exist for CRFs [START_REF] Szeliski | A comparative study of energy minimization methods for Markov random fields[END_REF]. One of the fastest is the α-expansion algorithm of [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], which is relying on graph-cut algorithms of [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF].

In the sparsity literature a number of algorithms have been proposed to take advantage computationally of the sparsity of the solution. In the convex setting notable examples are the LARS homotopy algorithms [START_REF] Efron | Least angle regression[END_REF] to compute the regularization path of the Lasso problem or working set algorithms such Glmnet [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF] that computes efficiently the solutions to 1 regularized problems. The Frank-Wolfe algorithm has also be used to exploit the sparsity of solution of optimization problem in the constrained setting [START_REF] Jaggi | Revisiting Frank-Wolfe: projectionfree sparse convex optimization[END_REF] as well as in the regularized setting [START_REF] Harchaoui | Conditional gradient algorithms for normregularized smooth convex optimization[END_REF].

In the non-convex setting greedy approaches that are used to compute sparse partial solutions count forward selection algorithms such as orthogonal matching pursuit [START_REF] Mallat | Adaptive timefrequency decomposition with matching pursuits[END_REF], orthogonal least squares [START_REF] Chen | Orthogonal least squares learning algorithm for radial basis function networks[END_REF] and related algorithms [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], and algorithms such as SBR [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] and FoBa [START_REF] Zhang | Adaptive forward-backward greedy algorithm for sparse learning with linear models[END_REF], which use backwards steps to remove previously introduced variables that are no longer relevant. See [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF] for a review.

A working set algorithm for total variation regularization

We propose to solve the minimization of a differentiable function f regularized by a weighted total variation of the form TV(x) = 1 2 ij∈E w ij |x i -x j |. Our working set algorithm alternates between solving a reduced problem of the form min x∈span(Π) Q(x) for Q(x) = f (x) + λTV(x), and refining the partition Π. We will discuss in Section 2.3 how to solve the reduced problem efficiently, but first present a criterion to refine the partition Π. The propositions of the next section are proved in the appendix.

Steepest binary cut

To obtain a new partition, and given the current solution x Π = arg min x∈span(Π) Q(x) of the reduced problem we consider updates of x of the form

x Π +h u B with u B = γ B 1 B -γ B c 1 B c
for some set B ⊂ V and some scalars h, γ B and γ B c such that u B 2 = 1. We discuss later (in section 2.2) how the choice of B leads to a new partition and focus first on a rationale to choose B. A natural criterion is to choose the set B such that

Q decreases the most in the direction of u B . If we let Q (x, v) = lim h→0 h -1 (Q(x+hv)-Q(x)) so that, when d ∈ R n is a unit vector, Q (x, d) denotes the directional derivative of Q at x ∈ R n in the direction d, then this criterion requires to solve min B⊂V Q (x Π , u B ). Note that since Q (x, 1 ∅ ) = 0 then min B⊂V Q (x, 1 B ) ≤ 0.
To further characterize Q we decompose the objective function. For a current value of x, we denote by S(x) := {(i, j) ∈ E | x i = x j } the set of edges connecting nodes with different values, which is arguably the appropriate notion of support for our setting, and S c (x) := E\S(x). Setting S := S(x Π ) and given that the absolute value is differentiable everywhere except at zero, we can split Q into two parts, Q S and TV |S c , which are respectively differentiable and non-differentiable at

x Π Q S (x) = f (x) + λ 2 (i,j)∈S w ij |x i -x j | TV |S c (x) = λ 2 (i,j)∈S c w ij |x i -x j |
TV |S c is a weighted total variation on the graph G but with weights w S c such that [w S c ] i,j = w ij for (i, j) ∈ S c and 0 else. We extend the previous notations and define

w S c (A, B) = w S c (A × B) = w((A × B) ∩ S c ). Proposition 1. For x ∈ R n , if we set S = S(x) then Q (x, 1 B ) = ∇Q S (x), 1 B +λw S c (B, B c ). Moreover if ∇f (x), 1 B = 0 then Q (x, u B ) = (γ B + γ B c ) Q (x, 1 B ).
Considering the case of x = x Π then for S = S(x Π ) clearly ∇Q S (x Π ) is orthogonal to span(Π) and thus to 1. Therefore, by the previous lemma, finding the steepest descent direction of the form u B requires to solve min

B⊂V (γ B + γ B c ) Q (x Π , 1 B ).
To keep a formulation which remains amenable to efficient computations, we will assume that γ B + γ B c is constant or ignore this factor2 . This leads us to define a steepest binary cut as any cut (B Π , B c Π ) such that

B Π ∈ arg min B⊂V ∇Q S (x Π ), 1 B +λw S c (B, B c ).
If ∅ is a solution, we set B Π = ∅. As formulated, this problem can be interpreted as a minimum cut problem in a suitably defined flow graph. Indeed consider G f low = (V ∪ {s, t}, E f low ) where s and t are respectively a source node and a sink node, and where the edge set E f low and the associated nonzero (undirected) capacities c ∈ R |S c |+n are defined as follows:

E flow =      (s, i), ∀i ∈ ∇ + c si = ∇ i Q S (x) (i, t), ∀i ∈ ∇ - c it = -∇ i Q S (x) (i, j), ∀(i, j) ∈ S c c ij = λw ij , (3) with ∇ + = {i ∈ V | ∇ i Q S (x) > 0} and ∇ -= V \∇ + . Proposition 2. Let S = S(x) then (C, V f low \C) is a minimal cut in G f low if and only if C\{s}, and its complement in V are minimizers of B → Q (x, 1 B ).
We can now characterize the optimality of the partition Π with its steepest binary partition B Π .

Proposition 3. We have x = arg min z∈R n Q(z) if and only if min B⊂V Q (x, 1 B ) = 0 and Q (x, 1 V ) = 0.
This guarantees that if we fail to find a steep cut, the algorithm has reached its optimum.

Induced new partition in connected sets

For Π = (A 1 , • • • , A k )
we motivated the choice of B Π by the best addition of a term for the form hu B to x = k i=1 c i 1 Ai . At the next iteration, we could thus consider minimizing Q under the constraint that x ∈ span(1 A1 , . . . ,

1 A k , 1 B ) with B = B Π . But, on this set, the values x i1 , x i2 , x i3 and x i4 with i 1 ∈ A j ∩ B, i 2 ∈ A j ∩ B c , i 3 ∈ A j ∩ B and i 4 ∈ A j ∩ B c
are a priori coupled; also, if A j ∩ B has several connected components, i → x i must take the same value on these components. These constraints are unnecessarily restrictive. Indeed, since the right notion of sparsity of the model is arguably the set of edges that can be non zeros S new = S(x) ∪ (B × B c ), it makes sense to only constrain the variables to stay in the set X Snew = {x | S(x ) = S new }. If we define Π new as the collection of all connected components in G of all sets A j ∩ B Π and A j ∩ B c Π for A j ∈ Π then it can be easily shown that span(Π new ) = X Snew .

Proposition 4. If B Π = ∅, Q(x Πnew ) < Q(x Π ). Algorithm 1: Working set scheme Initialize Π ← {V }, x Π ∈ arg min z=c1 V ,c∈R Q(z) while min B⊂V Q (x Π , 1 B ) < 0 do Pick B Π ∈ arg min B⊂V Q (x Π , 1 B ) Π ← {B Π ∩ A} A∈Π ∪ {B c Π ∩ A} A∈Π Π ← connected components of elements of Π Pick x Π ∈ arg min z∈span(Π) Q(z) return (Π, x Π )
We thus obtain Algorithm 1, illustrated on Figure 1. Based on the previous propositions, this algorithm guarantees monotonic convergence to the optimum Π . In terms of complexity, at each iteration Π new has at least one more component than Π, so that the algorithm converges in at most n steps. We now discuss how to exploit the sparse structure of x Π to solve the reduced problem efficiently.

A reduced graph for the reduced problem

Let Π be a coarse3 partition of V into connected components. We argue that min z∈span(Π) Q(z) can be solved on a weighted graph whose nodes are associated with the elements of the partition Π. Indeed, consider the graph

G = (V, E) with V = Π and E = {(A, B) ∈ V 2 | ∃(i, j) ∈ (A × B) ∩ E}. For x ∈ span(Π) we can indeed express TV(x) simply. Proposition 5. For x = A∈Π c A 1 A and 2TV G (c) := (A,B)∈E w(A, B) |c A -c B | we have TV(x) = TV G (c).
Note that if TV is the total variation associated with the weighted graph G with weights (w ij ) (i,j)∈E then TV G is the total variation associated with the weighted graph G and the weights w(A, B) (A,B)∈E . Denoting f : c → f ( A∈Π c A 1 A ), the reduced problem is equivalent to solving min c∈R k f (c) + λTV G (c) on G. Since Π is coarse, we have |E| m and hence computations involving TV G are much cheaper than those involving TV. As discussed in Section 2.4 the structure of f can often be exploited as well. G is cheap to compute comparatively to the speed ups allowed, as it is obtained by computing the connected components of the graph (V, E\S(x)), which can be done in linear time by depth-first search.

Solving linear inverse problems with TV

A number of classical problems in image processing such as deblurring, blind deconvolution and inpainting are formulated as ill-posed linear inverse problems [START_REF] Chan | Recent developments in total variation image restoration[END_REF], where a low TV prior on the image provides appropriate regularization. Typically if x 0 is the original image, H a blurring linear operator typically computed as a convolution on the image, additive noise and y = Hx 0 + the degraded image, this leads to problems of the form

x = arg min x∈R n 1 2 Hx -y 2 + λTV(x). ( 4
)
Since n is large, manipulating directly the matrices H or H should be avoided, but if the forward operator x → Hx is a convolution, it can be computed quickly using e.g. the fast Fourier transform. In the case of a blurring operator with adequate symmetry H Hx is also a blurring operator. These fast computations of x → Hx can be exploited in our algorithm to compute very efficiently the Hessian and loading vector of the reduced problem. Indeed, For a k-partition Π of V we denote by K ∈ {0, 1} n,k the components matrix whose columns are the vectors 1 A for A ∈ Π. Any x ∈ span(Π) can be rewritten as Kc with c ∈ R k . The gradient of the discrepancy function with respect to c writes: ∇ c 1/2 HKc -y 2 = K H HKc -K Hy As a result, the reduced problem can be solved by a similar forward backward scheme of much smaller size, with parameter K H HK and K Hy, which are of size k × k and k respectively.

Minimal partition problems

We consider now a generalization of the minimal partition problem of the form min

x∈R n Q(x) with Q(x) = f (x) + λΓ(x) where Γ(x) = 1 2 (i,j)∈S(x) w ij
the Mumford-Shah penalty. This non-convex nondifferentiable problem being significantly harder than the previous one, we restrict the functions f we consider to be separable functions of the form f (x) = i∈V f i (x i ) with f i : R → R continuously differentiable and convex. Inspired by greedy algorithms in the sparsity literature such as OMP, and by the working set algorithm we presented for TV regularization, we propose to exploit that |Π * | is not too large to construct an algorithm that greedily optimizes the objective by adding and removing cuts in the graph.

A greedy algorithm

As for the working set algorithm, we propose to gradually build an expansion of x of the form x = k i=1 c i 1 Ai for Π = (A 1 , • • • , A k ) a partition of V , by recursively splitting some of the existing sets A ∈ Π. As a first remark, given the separability assumption on f the choice of a global binary cut of a current partition Π reduces to cut individually all elements of Π, and the optimal cut for each A ∈ Π considered does furthermore not depend on x A c or Π\{A}. We should thus focus on cutting a single set A at a time. Without loss of generality, we consider the case Π = {V }.

Optimal binary cut

The optimal binary partition (B,

B c ) of V is such that Q optimized over span(1 B , 1 B c ) is as small as possi- ble. Since Γ(h1 B + h 1 B c ) = Γ(1 B ) = w(B, B c ) the corresponding optimization problem is of the form min B⊂V, h,h ∈R i∈B f i (h) + i∈B c f i (h ) + λw(B, B c ). ( 5
)
This problem is a priori hard to solve in general

(B → min h,h ∈R f (h1 B + h 1 B c ) is not submodular).
However, when h, h are fixed the assumption that f is separable entails that B → f (h1 B + h 1 B c ) is a modular function, so that the objective can be optimized with respect to B by solving a max-flow problem similar to (3) where f i (x i ) is substituted for ∇ i Q S (x). The smoothness and convexity of f w.r.t. h and h guarantee that the objective can be minimized efficiently w.r.t. these variables. A local minimum of the objective can thus be obtained efficiently by alternatively minimizing with respect to B and (h, h ) as suggested by [START_REF] Bresson | Fast global minimization of the active contour/snake model[END_REF] or El-Zehiry et al. (2011). One issue is again the fact that the binary partition B obtained as a solution of ( 5) and its complement are not necessarily connected sets, and it would be intractable to optimize only over pairs of sets (B, B c ) that are connected. However, splitting B and B c into their connected components will not increase the contour length Γ and can only decrease f . Consequently, given the collection of connected components A 1 , . . . , A k of B and B c we set

x = h 1 1 A1 + . . . + h k 1 A k with h j the minimizer of h → i∈Aj f i (h).

Recursive splitting and merging

Like for the working set algorithm, we recursively split the components of the current partition Π, however, in this case, as discussed earlier, only one connected component is cut at a time and the new boundary introduced is within that component. Our algorithm is structured by the following observations: Optimal cuts of existing components are independents. Given the separability assumption on f the problem of the choice of an optimal cut in each component are independent and need to be recomputed only for the newly created components.

Merging components. We consider a form of backward step, which merges a pair (A, B) of elements of Π and therefore removes the boundary between A and B. The step is taken if the increase of f induced is smaller than λw(A, B). This step is similar to the backward step of the Single Best Replacement algorithm by [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] which considers at each iteration the addition or removal of a single variable, whichever reduces most the value of the objective. Our approach differs in the sense that the merge we consider has no reason to correspond to a previous cut, since we consider all pairs of adjacent components. Saturated components. If the optimal binary partition of a component is the empty set, this component is said to be saturated. The separability of the fidelity function ensures that this component will not longer be split and can be ignored for the rest of the computation, unless it is involved in a mergeing step. The resulting algorithm, 0 Cut Pursuit is detailed in the supplementary material.

Experiments

Deblurring experiments with TV

To assess the performance in terms of speed of our working set algorithm we compare it with several state-of-the-art algorithms on a deblurring task of the form presented in section 2.4. Specifically, given an image x, we compute y = Hx + , where H is a Gaussian blur matrix, and is some Gaussian additive noise, and we solve (1) with a total variation regularization based on the 8-neighborhood graph built on image pixels. We use three 512 × 512 images of increasing complexities to benchmark the algorithms: the Shepp-Logan phantom, a simulated example and Lena represented in Figure 2 and 5 (all figures are represented in full size in the appendix). For all images the standard deviation of the blur is set at 5 pixels.

Competing methods

Preconditioned Generalized Forward Backward (PGFB). As a general baseline, we consider a recent preconditioned generalized forward-backward splitting algorithm by [START_REF] Raguet | Preconditioning of a generalized forward-backward splitting and application to optimization on graphs[END_REF] whose prior non-preconditioned version was shown to outperform state-of-the art convex optimization on deblurring tasks in [START_REF] Raguet | A generalized forward-backward splitting[END_REF].

Accelerated FB with parametric max-flows (FB+). Since efficient algorithms that solve the ROF problem have been the focus of recent work, and given that the ROF problem corresponds to the computation of the proximal operator of the total variation, we also make a comparison with an implementation of the accelerated forward-backward algorithm of [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]. To compute the proximal operator, we use an efficient solver of the ROF problem based on a reformulation as a parametric max-flow proposed by [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF]. The solver we used is the one made publicly available by the authors.

Cut Pursuit with Frank-Wolfe descent directions (CP-FW). We consider an alternative to the steepest binary partition to split the existing components of the partial solution. Inspired by a version of the Frank-Wolfe algorithm for regularized problems proposed by [START_REF] Harchaoui | Conditional gradient algorithms for normregularized smooth convex optimization[END_REF], we add the cut corresponding to the point of the set {x ∈ R n | Φ(x) ≤ 1} with minimal dot product with the gradient direction, which is itself computed as the solution of a maximum flow problem.

Cut Pursuit. To implement our algorithm (CP), we solve the min-cut problems with the solver of [START_REF] Kohli | Efficiently solving dynamic Markov random fields using graph cuts[END_REF] which itself is based on [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] and [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. The problems on the reduced graph are solved using the PGFB algorithm. This last choice is motivated by the fact that the preconditioning is quite useful as it compensates for the fact that the weights on the reduced graph can be quite imbalanced.

Results

Figure 3 shows the speed of the different algorithms on the three test images on a quad core CPU at 2.4 Ghz. We represent the relative primal suboptimality gap (Q t -Q ∞ )/Q ∞ where Q ∞ is the lowest value obtained by CP in 100 seconds. We can see that our algorithm speeds up significantly the direct optimization approach (PGFB) when the solution is sparse, and that it remains competitive in the case of a natural image with strong regularization. Indeed since the reduced problems are of much smaller size than the original, our algorithm can perform many more forward-backward iterations in the same allotted time.

The variant of Cut Pursuit using Frank-Wolfe directions (CP-FW) is as efficient over the first few iterations but then stagnates. The computation of a new Frank-Wolfe direction does not take into account the current support S(x), that provides a set of edges that are "free"; this entails that it overestimates the cost of adding new boundaries, resulting in too conservative updates. Accelerated forward-backward with parametric max-flow FB+ is slower in this setting, which is explained by the small number of updates the algorithm can afford in the allotted time given the cost of solving parametric max-flows; PGFB can make more than 50 updates during each iteration of FB+.

We report and discuss in the appendix the breakdown of computation time for each algorithm. Most notably the actual time spent by Cut Pursuit solving the reduced problem (with PGFB) takes comparatively very little time (around 3%) when this is the only step actually decreasing the objective function.

Experiments on minimal partitions

Denoising experiment

We now present experiments demonstrating the efficiency of the 0 -Cut Pursuit presented in section 3. We assess its performance against two state-of-the-art algorithms to minimize the Mumford-Shah energy of two noisy 512 × 512 images: the phantom of [START_REF] Shepp | The Fourier reconstruction of a head section[END_REF] and another simulated example. In order to illustrate the advantage of our algorithm over alternatives that discretize the value range, we make a small random shift of grey values to both images.

We also test the algorithms on a spatial statistic aggregation problem which corresponds to the design of a simplified map of the population density in the Paris area. We use raster open-source data 4 producing the map of Figure 5. The raster is partitionned using constrained Delaunay triangulation [START_REF] Chew | Constrained Delaunay triangulations[END_REF] to obtain a graph with 252,183 nodes and 378,258 edges. We use the squared loss weighted by the surface of each triangle as a fidelity term.

Competing methods

α-expansions on quantized models. Following [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF], we use the fact that if the range of values of x i is discretized, the MPP and TV problems reduce to maximum a posterior inference in a Potts model. More precisely it is a Potts model in which each value taken by the class variable c i is associated with a (non necessarily connected) level-set, and that has pairwise terms of the form 1 {ci =cj } w ij . We use α-expansions [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] to minimize the corresponding energy. More precisely we use the α-expansions implementation of [START_REF] Fulkerson | Class segmentation and object localization with superpixel neighborhoods[END_REF], which uses the same max-flow code [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF] as our algorithm. We denote the resulting algorithm (CRFi) where i is the number of quantization levels.

Non-convex relaxation. We implemented a nonconvex relaxation of the Mumford-Shah functional which is a "concave" version of the total variation, such as the adaptive Lasso [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] with t → ( + t) in lieu of t → |t|. The resulting functional can be minimized locally using a reweighted TV [START_REF] Ochs | On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision[END_REF]. We use our Cut Pursuit algorithm to solve each reweighted TV problems, since it was the fastest implementation. This algorithm is denoted TV 0.5 .

Results

We report on Figure 4 the energy obtained by the different algorithms normalized by the energy of the best constant approximation. We can see that our algorithms find faster local optima that are essentially as good or better than the ones obtained using αexpansions on the discretized problem, as long as the solutions are sufficiently sparse. For the population density data our algorithm performs roughly as well as CRF40 but is outperformed by CRF60.

Conclusion

We proposed two algorithms to minimize functions penalized respectively by the total variation and by the Mumford-Shah boundary size. They exploit that the solution typically has a sparse gradients, which entails that it has a small number of connected level-sets. constructing a sequence of approximate solutions that have the same property, they operate on reduced problems that can be solved efficiently, and require only to perform a number of graph cuts on the original graph, which are the bottleneck for further speed-ups. Like other working set schemes, our algorithms are not competitive if the solution has too many connected levelsets. In future work, we intend to extend the approach presented in this paper to approximate efficiently the total variation regularization path. It would also be interesting to find guarantees similar toS those existing for α-expansions to our greedy algorithm.

Figure 1 :

 1 Figure 1: Two first iterations of cut pursuit for the ROF problem on image (a). Image (b) and (d) represent the new cut at iterations 1 and 2, (c) and (e) the partial solution, with the current set of contours S in red.

Figure 2 :

 2 Figure 2: (a) Simulated example with Gaussian blur (b) deblurred image obtained by Cut Pursuit (CP).

Figure 5 :

 5 Figure 5: (a) Noisy Shepp-Logan phantom, (b) the result of 0 -CP, (c) triangulated population density in Paris area and (d) the simplified map obtained by 0 -CP (69% of the variance explained with 1.2% of the contour length).

4Figure 3 :

 3 Figure 3: Relative primal suboptimality gap Q t /Q ∞ -1 at time t (in seconds) for different algorithms on the deblurring task: FB+: , PGFB: , CP: CPFW: , for three 512 × 512 images and different regularization values: Shepp-Logan phantom (left), our simulated example (middle) and Lena (right).Each marker corresponds to a cut computation for CP, CPFW, and a proximal operator computation for FB+.

Figure 4 :

 4 Figure 4: Mumford-Shah energy at time t (in seconds) divided by the same energy for the best constant image obtained with: Non-convex relaxation (TV 0.5 ), 0 cut pursuit ( 0 -CP) and α-expansions (CRFi), for the Shepp-Logan phantom (left), our simulated example (middle) and the map simplification (right). Markers corresponds respectively to one reweighting, one α-expansion cycle and one cut for (TV 0.5 ), (CRF) and ( 0 -CP).

Figure 6 :

 6 Figure 6: PSNR at convergence and time to converge in seconds for the four algorithms as well as the noisy image for the first two denoising experiments.

Appearing in Proceedings of the 19 th International Conference on Artificial Intelligence and Statistics (AISTATS) 2016, Cadiz, Spain. JMLR: W&CP volume 41. Copyright 2016 by the authors.

Proximal methods that perform soft-thresholding or the non-convex IHT methods maintain sparse solutions as well, but typically need to update a full dimensional vector at each iteration, which is why we do not include them here. They blend however very well with active set algorithms.

γB and γBc could otherwise be determined by requiring that 1, uB = 0. More rigorously, descent directions considered could be required to be orthogonal to span(Π), but this leads to even less tractable formulations, that we therefore do not consider here.

A coarse partition corresponds to a sparse solution.


Timothée Jamin

Yacine Djama

Jean-Claude Bacri

Eric Falcon

Tuning the resonant frequencies of a drop by a magnetic field

Keywords: numbers: 47.65.Cb, 47.55.D-, 47.35.Pq

We report an experimental study of a magnetic liquid drop deposited on a superhydrophobic substrate and subjected to vertical vibrations in presence of a static magnetic field. It is wellknown that a flattened drop of usual liquid displays oscillating lobes at its periphery when vibrated. By adding ferromagnetic nanoparticles to a water drop and varying the strength of the magnetic field, we are experimentally able to efficiently tune the resonant frequencies of the drop. By using conservation energy arguments, we show that the magnetic field contribution is equivalent to adding an effective negative surface tension to the drop. Our model is found in good agreement with the experiments with no fitting parameter.

When a drop of liquid is dynamically driven by an external force, its free-surface generally displays an oscillating pattern at the drop resonant frequency. At a fundamental level, the dynamical study of such an oscillating drop occurs in various domains at different scales. Indeed, liquid drop behavior is used to model stellar mass collisions or neutron-star oscillations in astrophysics [START_REF] Durisen | Dynamic fission instabilities in rapidly rotating n = 3/2 polytropes -a comparison of results from finitedifference and smoothed particle hydrodynamics codes[END_REF], and nuclei in nuclear fission [START_REF] Amsden | Relativistic hydrodynamic theory of heavy-ion collisions[END_REF][START_REF] Zingman | Collective sow and pion production in a hydrodynamic model[END_REF], or in cellular biology [START_REF] Brangwynne | Active liquid-like behavior of nucleoli determines their size and shape in xenopus laevis oocytes[END_REF]. It is also involved in metrology for measuring viscosity of liquids [START_REF] Papoular | Gas-film levitated liquids: Shape fluctuations of viscous drops[END_REF] or surface tension of molten materials [START_REF] Fraser | Surface tension measurements on pure liquid iron and nickel by an oscillating drop technique[END_REF][START_REF] Przyborowski | Surface tension measurement of molten silicon by the oscillating drop method using electromagnetic levitation[END_REF]. At a more practical level, it includes printing (ink drop generation), mixing of drop clouds, droplet manipulations in microfluidic, optofluidic, and pharmaceutical industry [START_REF] Oh | Shaken not stirred -on internal flow patterns in oscillating sessile drops[END_REF]. Since the pioneer works of Rayleigh [START_REF] Rayleigh | On the Capillary Phenomena of Jets[END_REF], the resonant frequencies of a drop have been studied experimentally by means of different forcing mechanisms: water drop on a vibrating plate [START_REF] Yoshiyasu | Self-induced vibration of a water drop placed on an oscillating plate[END_REF][START_REF] Okada | Observation of the shape of a water drop on an oscillating Teflon plate[END_REF][START_REF] Noblin | Vibrations of sessile drops[END_REF][START_REF] Brunet | Star-drops formed by periodic excitation and on an air cushion -A short review[END_REF] or on a vibrating bath [START_REF] Dorbolo | Resonant and rolling droplet[END_REF], in microgravity [START_REF] Apfel | Free oscillations and surfactant studies of superdeformed drops in microgravity[END_REF], by acoustic or air-flow levitation [START_REF] Shen | Parametrically excited sectorial oscillation of liquid drops floating in ultrasound[END_REF][START_REF] Bouwhuis | Oscillating and star-shaped drops levitated by an airflow[END_REF], by levitating it on its own vapor (Leidenfrost drops) [START_REF] Holter | Vibrations of evaporating liquid drops[END_REF][START_REF] Snezhko | Pulsating gliding transition in the dynamics of levitating liquid nitrogen droplets[END_REF][START_REF] Piroird | Reshaping and capturing Leidenfrost drops with a magnet[END_REF], or by Lorentz force for a metal liquid drop [START_REF] Fautrelle | Free-surface horizontal waves generated by low-frequency alternating magnetic fields[END_REF][START_REF] Kocourek | Stability of liquid metal drops affected by a highfrequency magnetic field[END_REF]. A new challenge would be to accurately control and tune the free oscillations of a drop in a non-intrusive way. Notably, it would be of primary interest to be able to shift the resonant frequencies of the drop to avoid, for instance, some annoying frequency bands in practical situations.

In this Letter, we study the dynamics of a flattened drop of magnetic fluid deposited on a superhydrophobic substrate vertically vibrated in presence of a weak magnetic field of tunable amplitude. The fluid used is a ferrofluid constituted of a stable suspension of nanometric magnetic particles diluted in a carrier liquid (water). Above a critical acceleration of vibrations, the drop undergoes a parametric instability leading to an azimuthal pattern around the drop that oscillates at half the forcing frequency. A star-shaped drop is then observed made of several oscillating lobes. We observe that the resonant frequencies of the drop depend on the strength of the field. We show quantitatively and theoretically that we are able to shift these frequencies by tuning the effec-tive surface tension of the drop through the field. Indeed the contribution of the magnetic field is equivalent to the one of a negative surface tension. Including nanometric magnetic particles within a drop appears thus as a first step to be able to control and tune the free oscillations of a drop, such a magnetic fluid being known to shift the onsets of usual hydrodynamics instabilities (such as Kelvin-Helmholtz, Rayleigh-Taylor or Rayleigh-Plateau ones) [START_REF] Rosensweig | Ferrohydrodynamics[END_REF][START_REF] Blums | Magnetic Fluids[END_REF]. Finally, note that ferrofluid drop manipulation with a magnet is also an important task in microfluidics or lab on chip technology [START_REF] Bormashenko | New investigations on ferrofluidics: Ferrofluidic marbles and magnetic-field-driven drops on su-perhydrophobic surfaces[END_REF][START_REF] Sterr | Rolling ferrofluid drop on the surface of a liquid[END_REF][START_REF] Nguyen | Magnetowetting and sliding motion of a sessile ferrofluid droplet in the presence of a permanent magnet[END_REF], and in dynamic self-assembly phenomena [START_REF] Timonen | Switchable static and dynamic self-assembly of magnetic droplets on superhyrophobic surfaces[END_REF].

The experimental set-up is shown in Fig. 1. A volume of ferrofluid (V = 1 mL) is put on a plate. Due to gravity, the drop is flattened and looks like a puddle (radius, R = 9.7 mm, and thickness, h = 3.4 mm -see insets of Fig. 1). The plate was coated with a spray [START_REF]Ultra Ever Dry c spray -See[END_REF], to obtain a drop/substrate contact angle measured of roughly 165 • . Such superhydrophobic substrate is crucial to minimize the pinning force of the drop at the contact line, and thus to obtain much better reproducible experiments [START_REF] Celestini | Vibration of submillimetersize supported droplets[END_REF][START_REF] Chang | Substrate constraint modifies the rayleigh spectrum of vibrating sessile drops[END_REF]. The plate is then subjected to vertical sinusoidal vibrations by means of an electromagnetic shaker. The frequency and amplitude of vibrations are in the range 0 ≤ f e ≤ 52 Hz and 0 ≤ A ≤ 3 mm, respectively. An accelerometer is fixed below the plate to measure the normalized acceleration of vibrations Γ = Aω 2 e /g with ω e = 2πf e and g = 9.81 m.s -2 the acceleration of gravity. One has 0 ≤ Γ ≤ 4. The ferrofluid puddle on the plate is placed between two horizontal coils (38 cm in mean diameter), generating a vertical static magnetic field B in the range 0 ≤ B ≤ 100 G, 99% homogeneous in the horizontal plane [START_REF] Browaeys | Surface waves in ferrofluids under vertical magnetic field[END_REF]. The drop oscillations are recorded by means of a high-speed camera (Phan-tomV10) located above the drop, with a 1000 frames per second sampling and a 1600×1200 pixels resolution. The ferrofluid used is an ionic aqueous suspension synthesized with 12.4% by volume of maghemite particles (Fe 2 O 3 ; 7 ± 0.3 nm in diameter) [START_REF] Massart | The ferrofluid synthesis has been performed by D. Talbot, Laboratory LI2C[END_REF]. Nanoparticle diameter and magnetic field values are small enough to avoid sedimen-Typeset by REVT E X tation or agglomeration due to vibrations, gravity and magnetic fields [START_REF] Rosensweig | Ferrohydrodynamics[END_REF]. Ferrofluid properties are: density, ρ = 1550 kg.m -3 , surface tension, γ = 43 ± 3 mN.m -1 , magnetic susceptibility, χ(B) ≈ χ(B = 0) = 1 [START_REF]For our range of values of B, the magnetic susceptibility of the ferrofluid χ(B) varies less than 3%, thus one fixes χ(B) = χ(B = 0)[END_REF], magnetic saturation, M sat = 36 × 10 3 A.m -1 and dynamic viscosity, 1.4 × 10 -3 N.s.m -2 . Its capillary length is then l c = γ/(ρg) = 1.7 mm.

To accurately measure azimuthal drop oscillations, we use large flattened drops (R ≫ l c , see inset of Fig. 1). The Bond number, quantifying the ratio between gravity and capillary forces acting on the drop, is Bo ≡ (R/l c ) 2 ≃ 33. The ratio between magnetic and capillary forces is quantified by the magnetic Bond number Bo m ≡ B 2 R/(µ 0 γ) with R the puddle radius, µ 0 = 4π × 10 -7 H.m -1 the magnetic permeability of the vacuum. For our range of B, one has Bo m ≤ 2. Moreover, the magnetic effect is much smaller than gravity one (Bo m /Bo ≤ 0.06) meaning that the strength of B is weak enough to not deform the flattened region of the puddle at rest.

We first carry out experiments with no magnetic field (B = 0). Figure 2 shows the parametric instability of a ferrofluid puddle subjected to vertical sinusoidal vibrations. Above a critical acceleration of vibration Γ c , an azimuthal pattern is observed in the horizontal plane at the drop periphery: lobes oscillate radially at half the forcing frequency, f e /2. When f e is increased, the number n of oscillating lobes (mode number) increases from n = 2 to 9 as shown in Fig. 2

and movie (see Supplemental Material [35]).

The resonant frequencies f n of such an inviscid drop are independent of the nature of the forcing and arise from an interplay between inertia and surface tension effects. In the limit 2R ≫ h, the flattened drop shape (see inset of

d 2 a n (t) dt 2 + ω 2 n a n (t) = 0, ω 2 n = γ ρR 3 n(n 2 -1) (1)
with ω n = 2πf n . Now, applying a vertical sinusoidal vibration of the substrate forces parametrically the drop by modulating gravity g → g(t) = g [1 + Γ cos(ω e t)] with Γ = Aω 2 e /g. Besides, the thickness h of the puddle is usually given [START_REF] Yoshiyasu | Self-induced vibration of a water drop placed on an oscillating plate[END_REF][START_REF] Brunet | Star-drops formed by periodic excitation and on an air cushion -A short review[END_REF] by the quasi-static balance between gravity and capillary energies per unit volume (2γ/h ≃ ρgh/2), i.e. h ≃ 2l c = 2 γ/(ρg). The hypothesis of constant volume of the puddle V = πR 2 h then gives its radius R = V /(πh). Thus, a temporal modulation of gravity g(t) induces those of the puddle thickness h(t) ∼ g(t) -1/2 , radius R(t) ∼ h(t) -1/2 ∼ g(t) 1/4 , and eigenfrequencies ω 2 n (t) ∼ R(t) -3 ∼ g(t) -3/4 using Eq. ( 1). For weak acceleration of vibrations Γ ≪ 1, the parametric oscillator equation governing the amplitude of the lobes a n (t) then reads [START_REF] Yoshiyasu | Self-induced vibration of a water drop placed on an oscillating plate[END_REF]

d 2 a n (t) dt 2 + ω 2 n 1 - 3 4 Γ cos(ω e t) a n (t) = 0 (2)
Equation ( 2) is the Mathieu's equation whose solutions are marginality curves (or instability "tongues") separating stable zones (no deformation of the drop) and unstable zones where azimuthal standing waves at the drop periphery oscillates at half the forcing frequency (f e /2), near the resonant frequencies f n [START_REF] Mathews | Mathematical Methods of Physics[END_REF]. This corresponds to the parametric instability shown in Fig. 2.

These tongues of instability are displayed in Fig. 3 lobes at critical acceleration Γ c . Γ min c is the minimum critical acceleration of each tongue. Equation ( 2) predicts that the minimum of each marginality curve occurs at twice the eigenmode f e = 2f n . Notice that experimental values for f n are smaller than values given by Eq. ( 1). This difference (< 18%) may be explained by the fact that the quasi-static approximation yielding Eq. ( 2) is valid only for low-frequency excitation. At higher frequency, the existence of axisymmetric modes shifts the minima of the marginality curves [START_REF] Noblin | Triplon Modes of Puddles[END_REF][START_REF] Noblin | Vibrated sessile drops: Transition between pinned and mobile contact line oscillations[END_REF]. Indeed, we checked that this difference becomes negligible when we remove these axisymmetric modes by performing a control experiment with a superhydrophobic steady plate in contact with the drop top [40]. One can now wonder what is the effect of an applied vertical magnetic field B on the drop dynamics.

When B is increased for a fixed n, we observe that the instability tongue is shifted towards a lower frequency (see arrows in Fig. 3). The minimum of this curve and thus the eigenfrequency f n (B) are found to decrease with B. For n = 5, a relative shift of f n of 16% is observed between extreme values of B used. One defines the absolute shift of the eigenfrequency as ∆Ω n (B) ≡ ω 2 n (B = 0) -ω 2 n (B), taking thus positive values. We plot in the inset of Fig. 4, the frequency shift ∆Ω n (B) as a function of B for different n. We find that ∆Ω n (B) ∼ B 2 for our range of B regardless of n. All data in the inset of Fig. 4 are found to collapse on a single curve when plotting ∆Ω n (B) as a function of B 2 n 3 (not shown). In order to compare with the model described below, ∆Ω n (B) is then displayed in the main Fig. 4 as a function of

B 2 n(n 2 -1), noting that n(n 2 -1) ≈ n 3 for n ≥ 3.
This frequency shift is not due to a geometrical effect mediated by B such as a drop lengthening along the field direction [START_REF] Blums | Magnetic Fluids[END_REF][START_REF] Arkhipenko | Study of a magnetized fluid drops shape in a homogeneous magnetic field[END_REF][START_REF] Bacri | Instability of ferrofluid magnetic drops under magnetic field[END_REF]. Indeed, using Eq. ( 1), a decrease of the puddle radius R with B (up to 6% here) would lead to an increase of the resonant frequency, a situation opposite to our observations (see Fig. 3). Moreover, the results of Fig. 4 are found again when adding a steady plate above the drop. Finally, note that no clear dependence of the critical acceleration Γ min c on B is observed. An elegant way to understand the physical origin of the drop eigenfrequency shift with B is to balance energies involved in this system. In a first step, we assume B = 0 and follow Rayleigh's model [START_REF] Rayleigh | On the Capillary Phenomena of Jets[END_REF]. The flattened drop shape is approximated by a cylinder of fluid (see Fig. 5a). We denote S ⊥ its surface area normal to gravity and magnetic field (i.e. top plus bottom areas) and S the peripheral surface area. At rest, they are denoted S ⊥,0 and S ,0 . We consider small radial deformations of the peripheral surface of amplitude a n (t) ≪ R around an instantaneous radius R(t). In polar coordinates, this reads r(θ, t) = R(t) + a n (t) cos(nθ) (see Fig. 5b). h is assumed constant with time and thus also S ⊥ due to volume conservation. It is known, since Rayleigh [START_REF] Rayleigh | On the Capillary Phenomena of Jets[END_REF], that the radial deformation a n induces an increase in S (see Figs. 5a-b) as Thus, capillary energy increases by

∆S (t) ≡ S (t) -S ,0 = πha 2 n (t)(n 2 -1) 2R (3)
∆E c (t) = γ∆S (t). (4)
Besides, kinetic energy of the liquid reads

E k (t) = da n dt 2 πρhR 2 2n . (5)
The conservation of energy d(E k + E c )/dt = 0 then leads to Eq. ( 1) [START_REF] Rayleigh | On the Capillary Phenomena of Jets[END_REF].

Let us now introduce the magnetic energy E m in the Rayleigh's model. For a linearly permeable ferrofluid of volume V [START_REF]For our range of values of B, the magnetic susceptibility of the ferrofluid χ(B) varies less than 3%, thus one fixes χ(B) = χ(B = 0)[END_REF], one has E m = -B. V MdV /2, with B the external magnetic field and M the local ferrofluid magnetization [START_REF] Rosensweig | Ferrohydrodynamics[END_REF]. The determination of M needs to take into account the ferrofluid boundary conditions. For instance, for a plane parallel to B, one has E m = -χV B 2 /(2µ 0 ) (see Fig. 5c) and for a plane normal to B, E m = -χV B 2 /[2µ 0 (1 + χ)] (see Fig. 5d). Thus, the magnetic energy of a ferrofluid layer is smaller when B is parallel rather than normal to its surface.

For an arbitrary ferrofluid shape, M is nonuniform and an effective demagnetizing factor D is usually defined, with 0 ≤ D ≤ 1 (depending on the shape), such that [START_REF] Rosensweig | Ferrohydrodynamics[END_REF]

E m = - χV B 2 2µ 0 (1 + χD) . (6)
For volumes bounded by surfaces either parallel or normal to B, we can define the ratio between surface area normal to B and total surface area, r S ≡ S ⊥ /(S ⊥ + S ).

Using known theoretical values of D for different geometries [START_REF] Arrott | Micromagnetics of curling configurations in magnetically soft cylinders[END_REF][START_REF] Aharoni | Demagnetizing factors for rectangular ferromagnetic prisms[END_REF], we show in Fig. 6 that D ≈ r S over the whole range of aspect ratios. This means that E m decreases when the aspect ratio favors surfaces parallel to B, i.e. D → 0 when r s → 0. Black line is computed from model in [START_REF] Arrott | Micromagnetics of curling configurations in magnetically soft cylinders[END_REF] for a cylinder of axis aligned with B; red line (light gray) is computed from [START_REF] Aharoni | Demagnetizing factors for rectangular ferromagnetic prisms[END_REF] for a square rod aligned with B; is for a cube (D = rS = 1/3 [START_REF] Aharoni | Demagnetizing factors for rectangular ferromagnetic prisms[END_REF]); is for an infinitely elongated square rod normal to B (D = rS = 1/2 [START_REF] Aharoni | Demagnetizing factors for rectangular ferromagnetic prisms[END_REF]). Dashed line represents D = rS (slope 1).

We can then replace the demagnetizing factor D by r S in Eq. ( 6). The variation of S due to the presence of peripheral lobes induces a variation of E m through r S . Noticing that S ⊥,0 = 2πR 2 and S ,0 = 2πRh at rest, a first-order Taylor expansion in ∆S for small deformations (a n ≪ R) leads to ∆E m (t) = γ m ∆S (t) [START_REF] Przyborowski | Surface tension measurement of molten silicon by the oscillating drop method using electromagnetic levitation[END_REF] with

γ m = - χ 2 hB 2 4µ 0 (1 + χ + h/R) 2 , (8)
a quantity always negative. Then, using ∆S from Eq.

(3) and energies from Eqs (4), ( 5) and ( 7), the conservation of energy d(E k + E c + E m )/dt = 0 finally leads to the resonant frequencies of the ferrofluid drop

ω 2 n (B) = γ + γ m (B 2 ) ρR 3 n(n 2 -1) • (9)
Using Eq. ( 1) then leads to

ω 2 n (B = 0) -ω 2 n (B) = - γ m (B 2 ) ρR 3 n(n 2 -1) • (10)
The B 2 n(n 2 -1) scaling is in good agreement with the one found experimentally (see Fig. 4) as well as for the theoretical prefactor without fitting parameter. Notice that Eq. ( 9) includes the usual capillary contribution (γ) and a magnetic one (γ m ) that depends on B. The magnetic term thus plays the role of a negative surface tension (γ m < 0) that thus reduces the drop resonant frequencies. The magnetic field can be then used to tune the effective surface tension, γ ef f ≡ γ + γ m , and thus ω n . For our ranges of B, using ferrofluid properties and geometry, one has γ m ∈ [-8.5, 0] mN.m -1 , that is up to 20% of γ. The analogy with surface tension arises from ∆E m = γ m (B 2 )∆S . This means that an increase of the drop surface area parallel to B favors its magnetization and thus decreases E m since γ m < 0. Consequently, B has a stabilizing effect on the lobes. Finally, note that a model of the dynamics of a ferrofluid drop confined between two plates [START_REF] Langer | Dynamics of labyrinthine pattern formation in magnetic fluids[END_REF][START_REF] Jackson | Hydrodynamics of fingering instabilities in dipolar fluids[END_REF] mentioned such a possible negative surface tension effect, but requires χ ≪ 1 and thus cannot apply here where χ = 1.

To conclude, we have studied the dynamics of parametric oscillations of a centimetric ferrofluid drop on a superhydrophobic plate subjected to vertical sinusoidal vibrations and a constant magnetic field. By adding ferromagnetic nanoparticles to a water drop, we are able to shift significantly its eigenfrequencies by tuning the magnetic field strength. Using energy conservation, we extend the Rayleigh's model and show that the resonant frequency shift is well captured by our model with no fitting parameter. We also show that the magnetic field acts as a negative surface tension, and is a way to tune the effective surface tension of the drop. Finally, the weakness of the field strength and the small size of ferromagnetic particles are favorable to miniaturization to plan to control the oscillations of centimeter-to-micro-scale drop in a new non-intrusive way for potential applications.

FIG. 1 .

 1 FIG. 1. (Color online) Experimental setup. Bottom inset shows a schematic of the puddle at two phases of the vibration. Top inset: top view of the ferrofluid puddle at rest.

FIG. 2 .

 2 FIG. 2. Top view of the azimuthal pattern displayed around the drop as a function of fe. Above a critical acceleration of vibration, an azimuthal pattern is observed normal to the vibration direction. When the forcing frequency fe is increased (0 ≤ fe ≤ 52 Hz), the number n of lobes, oscillating radially at fe/2, increases from n = 2 to 9 (from left to right and top to bottom). See [35] for a movie. Magnetic field B = 0. Volume V = 1 mL

  in an acceleration/frequency phase space (Γ c /Γ min c vs f e ) for B = 0 (open-symbols). Data were obtained experimentally for various f e by increasing Γ until we observe

ForcingFIG. 3 .

 3 FIG. 3. (Color online) Phase diagram of normalized rescaled critical acceleration Γc/Γ min c vs. fe. Curves are marginality curves separating stable and unstable zones for different mode numbers n = 3 (♦), 4 (△) and 5 (•). No instability occurs for Γ ≤ Γc, whereas azimuthal pattern around the puddle occurs within tongues. Open symbols corresponds to B = 0, full symbols to B = 0: B = 44 (red, dark gray), 77 (black) and 99 G (green, light gray). When B is increased (see arrows), the tongue are shifted towards a lower frequency, for each mode n.

3 FIG. 4 .

 34 FIG. 4. (Color online) Eigenfrequency shift ω 2 n (B = 0)ω 2 n (B) as a function of B 2 n(n 2 -1) for different modes n =3 ( ), 4 ( ), 5 (•) and 7 ( ). Solid line is the prediction from the model of Eq. (10) with no fitting parameter. Inset: Unrescaled frequency shift vs. B. Dashed lines have a slope 2.

FIG. 5 .

 5 FIG. 5. (Color online) Schematic view of the puddle: (a) at rest, the puddle is considered as a cylinder of liquid, (b) when lobes appear, peripheral surface S increases as well as capillary energy. (c) and (d) illustrate magnetization within a ferrofluid film parallel or normal to B.

1 BFIG. 6 .

 16 FIG.6. (Color online) Theoretical demagnetizing factors D vs. surface ratio rS. Black line is computed from model in[START_REF] Arrott | Micromagnetics of curling configurations in magnetically soft cylinders[END_REF] for a cylinder of axis aligned with B; red line (light gray) is computed from[START_REF] Aharoni | Demagnetizing factors for rectangular ferromagnetic prisms[END_REF] for a square rod aligned with B; is for a cube (D = rS = 1/3[START_REF] Aharoni | Demagnetizing factors for rectangular ferromagnetic prisms[END_REF]); is for an infinitely elongated square rod normal to B (D = rS = 1/2[START_REF] Aharoni | Demagnetizing factors for rectangular ferromagnetic prisms[END_REF]). Dashed line represents D = rS (slope 1).

We thank D. Talbot for the ferrofluid synthesis, M. Berhanu, P. Brunet, M. Costalonga, and C. Laroche for fruitful discussions, and A. Lantheaume, Y. Le Goas and M.-A. Guedeau-Boudeville for technical help. T. J. was supported by the DGA-CNRS Ph.D. program. This work was partially financed by ANR Turbulon 12-BS04-0005.