Sylvain Bontemps

V Könyves

F Motte

N Schneider

P Didelon

V Minier

P Saraceno

D Ward-Thompson

Ph André

A Men'shchikov

J Di Francesco

G White

S Molinari

L Testi

A Abergel

M Griffin

Th Henning

P Royer

B Merín

R Vavrek

M Attard

D Arzoumanian

C D Wilson

P Ade

H Aussel

J.-P Baluteau

M Benedettini

J.-Ph Bernard

J A D L Blommaert

L Cambrésy

P Cox

A Di Giorgio

P Hargrave

M Hennemann

M Huang

J Kirk

O Krause

R Launhardt

S Leeks

J Le Pennec

J Z Li

P G Martin

A Maury

G Olofsson

A Omont

N Peretto

S Pezzuto

T Prusti

H Roussel

D Russeil

M Sauvage

B Sibthorpe

A Sicilia-Aguilar

L Spinoglio

C Waelkens

A Woodcraft

A Zavagno

From filamentary clouds

Keywords: formation -stars, luminosity function, mass function -stars, protostars -ISM, clouds -ISM, structure -submillimeter, ISM

published or not. The documents may come

Introduction

The Herschel Space Observatory (Pilbratt et al. 2010) offers a unique opportunity to improve our global understanding of the earliest phases of star formation. Here, we present first highlights from the Gould belt survey, one of the largest key projects with Herschel (cf. [START_REF] André | The Dusty and Molecular Universe: A Prelude to Herschel and ALMA[END_REF], based on extensive farinfrared and submillimeter mapping of nearby molecular clouds with both the SPIRE [START_REF] Griffin | From Darkness to Light[END_REF] and PACS (Poglitsch et al. 2010) bolometer cameras. This SPIRE/PACS imaging survey will cover the bulk of the nearest (d ≤ 0.5 kpc) cloud complexes in the Galaxy, which are mostly located in the Gould belt, a giant (∼ 700 pc × 1000 pc), flat structure inclined by ∼ 20 • to the Galactic plane (e.g., Guillout 2001). Since the ∼ 15 ′′ angular resolution of Herschel around λ ∼ 200 µm is adequate for probing individual (∼ 0.01-0.1 pc) star-forming cores up to ∼ 0.5 kpc away, the cloud complexes of the Gould belt correspond to the volume of Galactic space where Herschel imaging can be best used to characterize in detail the earliest stages of star formation.

The immediate observational objective of the Gould belt survey is to obtain complete samples of nearby prestellar cores and Class 0 protostars with well characterized luminosities, tempera-Send offprint requests to: pandre@cea.fr ⋆ Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. ⋆⋆ Appendix A and Figs. 3, 4 are only available in electronic form via http://www.edpsciences.org tures, and density profiles, as well as robust core mass functions and protostar luminosity functions, in a variety of star-forming environments. An order of magnitude more cold prestellar cores than already identified from the ground are expected to be found in the entire survey, which should allow us to derive an accurate prestellar core mass function (CMF) from the pre-browndwarf to the intermediate-mass range. Thanks to its high sensitivity and large spatial dynamic range, this Herschel survey can also probe, for the first time, the link between low-density cirrus-like structures in the interstellar medium (ISM) and compact self-gravitating cores. The main scientific goal is to elucidate the physical mechanisms responsible for the formation of prestellar cores out of the diffuse ISM, which is crucial for understanding the origin of the stellar initial mass function (IMF).

Our first results, obtained toward the Aquila Rift and Polaris clouds, are very promising (e.g., Könyves et al. 2010, Bontemps et al. 2010, Men'shchikov et al. 2010, Ward-Thompson et al. 2010, and Miville-Deschênes et al. 2010). As discussed in Sect. 4 of this paper, they suggest that prestellar cores result from the gravitational fragmentation of filaments in the cold ISM.

Herschel observations

The Herschel survey was designed to cover the densest portions of the Gould belt with SPIRE at 250-500 µm and PACS at 100-160 µm. The observational goal is to make a complete, homogeneous mapping of the A V > 3 regions with SPIRE and the A V > 6 regions with PACS, and representative areas at A V ∼ 1-3 levels with both instruments. The survey sensitivity (better than A V ∼ 1 at the 5σ level) will allow us to probe the structure of nearby molecular clouds down to the level of the interface with their atomic gas envelopes. The 15 target clouds span a range of physical conditions, from active, cluster-forming complexes to quiescent regions with lower or no star formation activity 1 .

Our Herschel mapping consists of two steps:

(1) A wide-field SPIRE/PACS survey of a total surface area ∼ 160 deg 2 using the so-called "parallel mode" with a scanning speed of 60 ′′ /s and the PACS 70 µm and 160 µm bands. In this first step, the main goal is to acquire adequate SPIRE 250-500 µm data for all of the target regions. The PACS data acquired simultaneously yield very useful information at 70 µm and 160 µm through most of the SPIRE survey, but do not have optimal (diffraction-limited) angular resolution. We selected the 70 µm filter for the blue band of PACS to obtain a good diagnostic of the presence of embedded protostars (cf. [START_REF] Dunham | [END_REF] throughout the mapped regions.

(2) A dedicated PACS-only survey of a total surface area ∼ 65 deg 2 , observing the 100 µm and 160 µm bands with a scanning speed of 20 ′′ /s. This second step is supposed to provide data of optimal sensitivity and resolution at 100 µm and 160 µm.

The Aquila Rift and Polaris flare regions were observed in parallel mode with both SPIRE and PACS (step 1) during the science demonstration phase of Herschel. The corresponding observations are described in detail by Könyves et al. (2010), Bontemps et al. (2010), Men'shchikov et al. (2010), and Ward-Thompson et al. (2010) in this special A&A issue. The Polaris flare field is a high-latitude translucent cloud with little to no star formation at d ∼ 150 pc (e.g., Heithausen et al. 2002), and is expected to have the lowest level of background cloud emission and cirrus confusion noise in the entire Gould belt survey. At the other extreme, the Aquila field is a very active star-forming complex at d ∼ 260 pc (e.g., [START_REF] Gutermuth | [END_REF], and is expected to have the highest level of background cloud emission and cirrus confusion noise in the whole survey. Sensitive Galactic far-IR imaging surveys such as the one discussed here are not limited by instrumental sensitivity but by confusion arising from small-scale cirrus/cloud structure (cf. Roy et al. 2010 and references therein). The levels of background cloud fluctuations observed in the Aquila and Polaris fields allow us to roughly estimate the whole range of cirrus confusion noise levels that will affect the Gould belt survey. In the Aquila field, the observed rms level of cirrus or "structure" noise (as measured on the typical scales of dense cores) ranges from ∼ 30 mJy/18 ′′beam to ∼ 300 mJy/18 ′′ -beam at λ = 250 µm, which is a factor of ∼ 3-30 above the rms instrumental noise. In contrast, in the Polaris field, the rms level of emission fluctuations measured in the SPIRE 250 µm map is only ∼ 10-30 mJy/18 ′′ -beam, which is very close to the rms sensitivity level expected after two crossscans in parallel mode at 60 ′′ /s. Assuming the dust opacity of Hildebrand (1983) and the median dust temperatures derived from our Herschel maps, the above surface-brightness sensitivity levels translate into 5σ column-density detection thresholds of N H2 < ∼ 3 × 10 20 cm -2 in Polaris and N H2 < ∼ 10 21 cm -2 in Aquila (in agreement with the column-density maps shown by Ward-Thompson et al. 2010 andKönyves et al. 2010). This high column-density sensitivity, coupled with the unprecedented surface brightness and spatial dynamic ranges of our Herschel mapping, allows us to probe, for the first time, the physical connection between the structure of the diffuse ISM and the formation of prestellar cores (cf. Sect. 4).

1 See http://gouldbelt-herschel.cea.fr/ for the list of target fields.

The corresponding mass sensitivities for typical prestellar/protostellar cores are ∼ 0.3 M ⊙ (85% completeness level) in Aquila (d = 260 pc) and ∼ 0.01 M ⊙ (85% completeness level) in Polaris (d = 150 pc). These completeness numbers were estimated by performing Monte-Carlo simulations as described by Könyves et al. (2010) in the case of Aquila. The initial results obtained in Aquila and Polaris therefore confirm that the completeness level of our Herschel census for prestellar cores will reach the pre-brown dwarf mass regime (M core < 0.08 M ⊙ ) in the nearest molecular clouds (d ∼ 150 pc) of the Gould belt.

Main results and analysis

The Herschel images of the Aquila Rift and Polaris flare regions exhibit extensive filamentary structure, as well as numerous dense cores situated along these filaments (see Fig. 1, online Fig. 3, and Men'shchikov et al. 2010). A total of 541 starless cores (∼ 0.01-0.1 pc in size) can be identified in the whole ( > ∼ 3.3 • × 3.3 • ) Aquila field, most (> 60%) of which appear to be self-gravitating and prestellar in nature. The latter is inferred from: (1) a comparison of the core masses derived from the SPIRE/PACS spectral energy distributions with local values of the Jeans or Bonnor-Ebert (BE) mass also estimated from Herschel data (see Könyves et al. 2010); (2) the positions of the cores in a mass versus size diagram (cf. online Fig. 4), which are close to the loci expected for critical isothermal BE spheres at gas temperatures ∼ 7-20 K; and (3) the mean column densities of the cores, which exceed the background column densities by a median factor ∼ 1.5 as expected for critically self-gravitating BE spheres. The shape of the core mass function (CMF) derived for this sample of Aquila cores closely resembles the IMF (Fig. 2-left -see Könyves et al. 2010 for the CMF of the Aquila central region and Appendix A for details on the derivation of core masses). In contrast, the 302 starless cores identified in the cirrus-like Polaris cloud with the same clump-finding algorithm (getsources -see Men'shchikov et al. 2010) appear to be mostly unbound (cf. online Fig. 4) and their mass distribution does not follow the IMF, with a peak at an order of magnitude smaller mass (Fig. 2-right). Only 5 of the Polaris cores are reasonably close to being gravitationally bound and thus possibly prestellar in nature (see Ward-Thompson et al. 2010). Between 45 and 60 Class 0 protostars are revealed by Herschel in the Aquila field (Bontemps et al. 2010), while not a single protostar is detected in the Polaris region.

The Aquila filaments harboring embedded protostars and/or large concentrations of prestellar cores are all characterized by higher column densities (N H 2 > ∼ 10 22 cm -2 ), suggesting that they are gravitationally unstable (see Fig. 1-left). In contrast, both the Polaris filaments and the quiescent, non-star-forming filaments observed in Aquila have much lower column densities (up to a few 10 21 cm -2 ), suggesting they are gravitationally stable (see Fig. 1-right and Sect. 4 below).

Discussion and conclusions

The Herschel early results obtained toward the Aquila starforming complex confirm that the prestellar CMF resembles the stellar IMF, using data with already a factor of ∼ 2 to 9 better counting statistics than earlier (sub)-millimeter continuum or near-IR extinction surveys (e.g., Motte et al. 1998, Johnstone et al. 2000, Enoch et al. 2006, Stanke et al. 2006, Nutter & Ward-Thompson 2007, Alves et al. 2007). The close resemblance between the Aquila CMF and the IMF in both shape and 2010) are shown as green stars and blue triangles, respectively. Note the good correspondence between the spatial distribution of the bound cores/protostars and the regions where the filaments are unstable to gravitational collapse. mass scale suggests that, as a class, the self-gravitating prestellar cores identified in far-IR/submilimeter continuum imaging studies such as the present Herschel survey may form stars on a one-to-one basis, with a fixed and relatively high local efficiency, i.e., ǫ core ≡ M ⋆ /M core ∼ 20 -40% in Aquila. This is consistent with theoretical models according to which the stellar IMF is in large part determined by pre-collapse cloud fragmentation, prior to the protostellar accretion phase (cf. Larson 1985;Padoan & Nordlund 2002;Hennebelle & Chabrier 2008). There are several caveats to this simple picture (cf. discussion in [START_REF] André | Structure Formation in Astrophysics[END_REF], and detailed analysis of the data from the whole Gould belt survey will be required to fully characterize the CMF-IMF relationship and, e.g., investigate possible variations in the efficiency ǫ core with environment. It is nevertheless already clear that one of the keys to the problem of the origin of the IMF lies in a good understanding of the formation process of prestellar cores, even if additional processes, such as rotational subfragmentation of prestellar cores into binary/multiple systems (e.g., [START_REF] Bate | [END_REF], probably also play an important role.

Our Herschel initial results also provide key insight into the core formation issue. They support an emerging picture (see also Myers 2009) according to which complex networks of long, thin filaments form first within molecular clouds, possibly as a result of interstellar MHD turbulence, and then the densest filaments fragment into a number of prestellar cores via gravitational instability. That the formation of filaments in the diffuse ISM represents the first step toward core/star formation is suggested by the filaments already being omnipresent in a diffuse, non-starforming cloud such as Polaris (cf. Fig. 1-right, Men'shchikov et al. 2010, andMiville-Deschênes et al. 2010). The second step appears to be the gravitational fragmentation of a subset of the filaments into self-gravitating cores. This observationally-driven scenario can be placed on a stronger footing by comparing our Herschel results to existing models of filamentary cloud fragmentation. Inutsuka & Miyama (1992, 1997) showed that an unmagnetized isothermal filament is unstable to axisymmetric perturbations if the line mass or mass per unit length, M line , of the filament is larger than the critical value required for equilibrium, M unmag line,crit = 2 c 2 s /G, where c s is the isothermal sound speed. Remarkably, the critical line mass only depends on gas temperature (Ostriker 1964) and is modified by only a factor of order unity for filaments with realistic levels of magnetization (Fiege & Pudritz 2000). Figure 1 shows maps of the mass per unit length, expressed in approximate units of the critical line mass, for the filaments detected by Herschel in Aquila and Polaris. These maps were constructed from the column density maps derived from our SPIRE/PACS images (see Appendix A) by multiplying the local column density of the filaments by their measured typical width (FWHM ∼ 14000 AU in Aquila and ∼ 9000 AU in Polaris, before deconvolution -cf. Men'shchikov et al. 2010). In Fig. 1, the critical line mass corresponding to a gas temperature of 10 K (M 10K,unmag line,crit ≈ 15 M ⊙ /pc) was adopted throughout the fields, but nearly identical results are obtained if the dust temperature maps derived from the Herschel images (cf. Bontemps et al. 2010) are used instead. The results (cf. Fig. 1) show that most (> 60%) of the bound prestellar cores identified in Aquila are concentrated in supercritical filaments with M line > M unmag line,crit . Furthermore, virtually all supercritical fil-Fig. 2. Core mass functions (blue histograms with error bars) derived from our SPIRE/PACS observations of the Aquila (left) and Polaris (right) regions, which reveal of total of 541 candidate prestellar cores and 302 starless cores, respectively. A lognormal fit (red curve) and a power-law fit (black solid line) to the high-mass end of the Aquila CMF are superimposed in the left panel. The power-law fit has a slope of -1.5 ± 0.2 (compared to a Salpeter slope of -1.35 in this dN/dlogM format), while the lognormal fit peaks at ∼ 0.6 M ⊙ and has a standard deviation of ∼ 0.43 in log 10 M. The IMF of single stars (corrected for binaries -e.g., Kroupa 2001), the IMF of multiple systems (e.g., [START_REF] Chabrier | The Initial Mass Function 50 years later[END_REF], and the typical mass spectrum of CO clumps (e.g., Kramer et al. 1998) are also shown for comparison. Note the remarkable similarity between the Aquila CMF and the stellar IMF, suggesting a ∼ one-to-one correspondence between core mass and star/system mass with M ⋆sys = ǫ M core and ǫ ≈ 0.4 in Aquila.

aments harbor prestellar cores and/or embedded protostars, in agreement with the view that they are collapsing and actively forming stars at the present time. In particular, this is the case for the Aquila main filament, which has M line > 5 × M unmag line,crit and is associated with a very rich protocluster (Serpens South -see [START_REF] Gutermuth | [END_REF], Bontemps et al. 2010). In contrast, the subcritical filaments with M line < M unmag line,crit are generally devoid of prestellar cores and protostars, which is consistent with the view that they are gravitationally stable, hence neither collapsing nor forming stars. All of the Polaris filaments appear to be in the subcritical regime, the maximum observed value of the stability parameter being M line /M unmag line,crit ∼ 0.45, and it is unclear whether they will evolve into the unstable regime or not.

It is noteworthy that the critical line mass approximately corresponds to a critical column density N H 2 ,crit ∼ 10 22 cm -2 , i.e., to a critical visual extinction A V,crit ∼ 10. Our Herschel findings thus provide an explanation of the visual extinction threshold for core formation found by earlier ground-based studies (e.g., Onishi et al. 1998, Johnstone et al. 2004). Prestellar cores are only observed above a threshold A V,crit because they form out of a filamentary background and only the supercritical, gravitationally unstable filaments are able to fragment into bound cores.

Confirming and refining this scenario for core formation will require the results of the entire Herschel survey, as well as follow-up (sub-)millimeter dust polarimetry and molecular line observations to, e.g., clarify the roles of magnetic fields, turbulence, and gravity in forming the filaments. Our initial results are nevertheless extremely encouraging. Extrapolating from the > ∼ 500 prestellar cores and > ∼ 45 Class 0 protostars identified in the ∼ 19 deg 2 covered by the Aquila and Polaris fields, we expect that the ∼ 160 deg 2 Gould belt survey will reveal a total of about 4500 prestellar cores, including a large number of candidate pre-brown dwarfs in the nearest (d ∼ 150 pc) clouds, and more than 350 Class 0 protostars. This will provide a unique database, including the southern hemisphere, for follow-up high-resolution molecular line/dust continuum studies of the physics of individual cores and protostars with ALMA. . The masses were derived as explained in Appendix A (see also Könyves et al. 2010) and the sizes were measured at 250 µm. For reference, the locations of the (sub)mm continuum prestellar cores identified by Motte et al. (1998Motte et al. ( , 2001) ) in ρ Oph and NGC2068/2071, respectively, are shown, along with the correlation observed for diffuse CO clumps (shaded band -cf. Elmegreen & Falgarone 1996). The two black solid lines mark the loci of critically self-gravitating isothermal Bonnor-Ebert spheres at T = 7 K and T = 20 K, respectively. The 341 cores classified as prestellar by Könyves et al. (2010), out of a total of 541 starless cores in Aquila, are shown as filled triangles in the left panel. The 302 starless cores of Polaris lie much below the two Bonnor-Ebert lines in the right panel, suggesting they all are unbound. The red solid lines are two lines of constant mean column density at N H 2 = 10 21 cm -2 and N H 2 = 10 22 cm -2 , respectively. The 5σ detection threshold at d = 150 pc of existing ground-based (sub)mm (e.g., SCUBA) surveys as a function of size is shown by the dashed curve. (The shape of this curve reflects a constant sensitivity to column density until source size approaches the beam size.) The 5σ detection thresholds of our SPIRE 250 µm observations, given the estimated levels of cirrus noise (cf. Sect. 2) and assumed distances of 260 pc for Aquila and 150 pc for Polaris, are shown by the blue curves. The blue arrow in the left panel indicates how the positions of the Herschel cores (and the blue curve marking the detection threshold) would move in the diagram if a distance of 400 pc were adopted for Aquila instead of 260 pc (see Appendix A).

Fig. 1 .

 1 Fig. 1. Column density maps of two subfields in Aquila (left) and Polaris (right) derived from our SPIRE/PACS data. The contrast of the filaments with respect to the non-filamentary background has been enhanced using a curvelet transform as described in Appendix A. Given the typical width ∼ 10000 AU of the filaments, these column density maps are equivalent to maps of the mass per unit length along the filaments. The color scale shown on the right of each panel is given in approximate units of the critical line mass of Inutsuka & Miyama (1997) as discussed in Sect. 4. The areas where the filaments have a mass per unit length larger than half the critical value and are thus likely gravitationally unstable have been highlighted in white. The maximum line mass observed in the Polaris region is only ∼ 0.45 × the critical value, suggesting that the Polaris filaments are stable and unable to form stars at the present time. The candidate Class 0 protostars and bound prestellar cores identified in Aquila by Bontemps et al. (2010) and Könyves et al. (2010) are shown as green stars and blue triangles, respectively. Note the good correspondence between the spatial distribution of the bound cores/protostars and the regions where the filaments are unstable to gravitational collapse.

Fig. 3 .

 3 Fig. 3. Composite 3-color images of the ∼ 11 deg 2 Aquila field (left) and the ∼ 8 deg 2 Polaris field (right) produced from our PACS/SPIRE parallel-mode images at 70, 160, and 500 µm. The color coding of the Aquila image is such that Red = SPIRE 500 µm, Green = PACS 160 µm, Blue = PACS 70 µm. For the Polaris image, Red = combination of the three SPIRE bands, Green = PACS 160 µm, Blue = PACS 70 µm. The Aquila composite image was also the first release of "OSHI", ESA's Online Showcase of Herschel Images (cf. http://oshi.esa.int and http://www.esa.int/SPECIALS/Herschel/SEMT0T9K73G 0.html).

Fig. 4 .

 4 Fig.4. Mass vs. size diagrams for the starless cores detected with Herschel-SPIRE/PACS in Aquila (left) and Polaris (right) (blue triangles). The masses were derived as explained in Appendix A (see alsoKönyves et al. 2010) and the sizes were measured at 250 µm. For reference, the locations of the (sub)mm continuum prestellar cores identified byMotte et al. (1998Motte et al. ( , 2001) ) in ρ Oph and NGC2068/2071, respectively, are shown, along with the correlation observed for diffuse CO clumps (shaded band -cf.Elmegreen & Falgarone 1996). The two black solid lines mark the loci of critically self-gravitating isothermal Bonnor-Ebert spheres at T = 7 K and T = 20 K, respectively. The 341 cores classified as prestellar byKönyves et al. (2010), out of a total of 541 starless cores in Aquila, are shown as filled triangles in the left panel. The 302 starless cores of Polaris lie much below the two Bonnor-Ebert lines in the right panel, suggesting they all are unbound. The red solid lines are two lines of constant mean column density at N H 2 = 10 21 cm -2 and N H 2 = 10 22 cm -2 , respectively. The 5σ detection threshold at d = 150 pc of existing ground-based (sub)mm (e.g., SCUBA) surveys as a function of size is shown by the dashed curve. (The shape of this curve reflects a constant sensitivity to column density until source size approaches the beam size.) The 5σ detection thresholds of our SPIRE 250 µm observations, given the estimated levels of cirrus noise (cf. Sect. 2) and assumed distances of 260 pc for Aquila and 150 pc for Polaris, are shown by the blue curves. The blue arrow in the left panel indicates how the positions of the Herschel cores (and the blue curve marking the detection threshold) would move in the diagram if a distance of 400 pc were adopted for Aquila instead of 260 pc (see Appendix A).

Acknowledgements. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC (UK); and NASA (USA). PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain).

Appendix A: Derivation of core/filament properties and effects of distance uncertainties

As described in more detail in a companion paper by Könyves et al. (2010) on Aquila, the masses of the cores identified in the Herschel images with the getsources algorithm (see Men'shchikov et al. 2010) were derived by fitting greybody functions to the spectral energy distributions (SEDs) constructed from the integrated flux densities measured with SPIRE/PACS for each core. We assumed the dust opacity law κ ν = 0.1 (ν/1000 GHz) 2 cm 2 /g, where ν denotes frequency and κ ν is the dust opacity per unit (gas + dust) mass column density. This dust opacity law, which is very similar to that advocated by Hildebrand (1983), is consistent with the value κ 1.3mm = 0.005 cm 2 /g adopted for starless cores in numerous earlier studies (e.g., Motte et al. 1998). Ignoring any systematic distance effect (see below), the core mass uncertainties are dominated by the uncertainty in κ ν , typically a factor of ∼ 2. Cores were classified as either protostellar or starless based on the presence or absence of significant PACS emission at 70 µm, respectively (cf. Bontemps et al. 2010 and[START_REF] Dunham | [END_REF]). In the Polaris field, the cirrus noise level is so low (cf. Fig. 4-right) that we cannot exclude that a fraction of the 302 candidate starless cores extracted with getsources correspond to background galaxies.

A column density map was derived for each region from the Herschel images smoothed to the SPIRE 500 µm resolution (36.9 ′′ FWHM) using a similar SED fitting procedure on a pixel by pixel basis (see, e.g., Figs. 1 and 6 of Könyves et al. 2010 for Aquila). To obtain the maps of the filamentary background shown in Fig. 1 of this paper, we then performed a "morphological component analysis" decomposition (e.g., Starck et al. 2003) of the original column density maps on curvelets and wavelets. The curvelet component images shown in Fig. 1 provide a good measurement of the column density distribution of the filamentary background after subtraction of the compact sources/cores since the latter are contained in the wavelet component. We estimate that these column density maps are accurate to within a factor of ∼ 2. The scaling in terms of the mass per unit length along the filaments is more uncertain, however, as it depends on distance (see below) and would in principle require a detailed analysis of the radial profiles of the filaments, which is beyond the scope of the present letter. Here, we simply assumed that the filaments had a Gaussian radial column density profile and multiplied the surface density maps by 2π 8 ln2 × W ≈ 1.06 W, where W is the typical FWHM width of the filaments. We assumed a mean molecular weight of µ = 2.33. At this stage, the correspondence between the critical line mass of the filaments, M unmag line,crit , and the visual extinction threshold, A V,crit ∼ 10 (see Sect. 4), is thus accurate to at best a factor of > ∼ 2.

There is some ambiguity concerning the distance to the Aquila Rift region. A number of arguments, presented in a companion paper by Bontemps et al. (2010), suggest that the whole region corresponds to a coherent cloud complex at d -= 260 pc (see also [START_REF] Gutermuth | [END_REF], which is the default distance adopted in the present paper for Aquila. However, other studies in the literature (see references in Bontemps et al. 2010) place the complex at a larger distance, d + = 400 pc. It is thus worth discussing briefly how our Aquila results would be affected if we adopted the larger distance estimate, d + , instead of d -. The core mass estimates, which scale as S ν d 2 /[B ν (T d ) κ ν ] where S ν is integrated flux density and B ν (T d ) is the Planck function, would systematically increase by a factor of 2.4. This would shift the CMF shown in Fig. 2-left to the right and thus lower the efficiency ǫ core from ∼ 20-40% to ∼ 10-20%. In the mass versus size diagram of Fig. 4, the cores would move upward as indicated in the left panel of the figure, which would increase the fraction of candidate bound cores in Aquila from 63% to 81%. The column density map of the Aquila filaments shown in Fig. 1 would be unchanged, but the scaling in terms of the mass per unit length along the filaments would change by ∼ 50% upward, since the physical width of the filaments would increase by ∼ 50%. In other words, the highlighted regions in Fig. 1left, where the mass per unit length of the filaments exceeds half the critical value, would slightly expand, increasing the contrast with the Polaris filaments and improving the correspondence between the spatial distribution of the prestellar cores/protostars in Aquila and that of the gravitationally unstable filaments. To summarize, our main conclusions do not depend strongly on the adopted distance.


Valentina Sessa
email: valentina.sessa@uerj.br

Luigi Iannelli
email: luigi.iannelli@unisannio.it

Francesco Vasca
email: vasca@unisannio.it

Vincent Acary
email: vincent.acary@inria.fr

A complementarity approach for the computation of periodic oscillations in piecewise linear systems

Keywords: Piecewise linear systems, Complementarity systems, Periodic solutions, Autonomous systems

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Piecewise linear (PWL) models can represent a wide class of practical systems and can exhibit interesting nonlinear behaviours such as periodic steady-state oscillations. In this paper, we consider the problem of the computation of periodic solutions for linear timeinvariant dynamical systems whose inputs are related to the outputs through a static multi input-multi output PWL relation, possibly set-valued. That class of dynamical systems has attracted a considerable interest in the literature, and it allows to capture even complex behaviours arising in many practical systems. We can mention electrical circuits where the current-voltage characteristics of some devices can be considered as piecewise linear, such as nonlinear resistors [START_REF] Stern | Piecewise-linear network theory[END_REF], ideal diodes and switches [START_REF] Camlibel | Switched networks and complementarity[END_REF][START_REF] Acary | Time-stepping numerical simulation of switched circuits within the nonsmooth dynamical systems approach[END_REF], or mechanical systems with Coulomb friction [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]. In those applications, the steadystate behaviour is generally depicted by a periodic motion [START_REF] Melin | Two types of limit cycles of a resonant converter modelled by threedimensional system[END_REF][START_REF] Yang | Limit cycle and chaotic invariant sets in autonomous hybrid planar systems[END_REF]. Since they tend to periodically oscillate also when no external excitation is applied, the period related to the periodic solution is difficult to be predicted a priori.

A considerable literature deals with the problem of computing steady-state periodic solutions and the corresponding period. The most popular time-domain technique is the so-called shooting method [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF][START_REF] Flieller | Computation and stability of limit cycles in hybrid systems[END_REF][START_REF] Li | A new method to determine the periodic orbit of a nonlinear dynamic systems and its period[END_REF]. A drawback of this method is the cost of calculating the sensitivity matrix, which is a crucial issue [START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF][START_REF] Theodosiou | On periodic steady state response and stability of Filippov-type mechanical models[END_REF][START_REF] Bizzarri | Steady state computation and noise analysis of analog mixed signal circuits[END_REF]. For instance, in [START_REF] Leine | Stick-slip vibrations induced by alternate friction models[END_REF][START_REF] Van De Vrande | An approximate analysis of dry-friction-induced stick-slip vibrations by smoothing procedure[END_REF], in order to compute periodic stick-slip vibrations belonging to mechanical systems with discontinuous dry friction, the shooting method is applied only after the approximation of the discontinuous characteristic by a suitable smooth function. Instead, in this paper, we show that the proposed representation allows to include the discontinuity into the model. For hybrid systems that include such class, in [START_REF] Thota | TC-HAT ( T C): a novel toolbox for the continuation of periodic trajectories in hybrid dynamical systems[END_REF] it is proposed an algorithm for the bifurcation analysis of periodic solutions. In that paper, a multi-point boundary-value problem is formulated by using the knowledge of the sequence of modes. The specific sequence is exploited also for fixing an initial guess for the periodic solution to be computed. In our paper, the periodic steady-state solution is obtained without fixing a priori the sequence of modes and without this information for choosing an initial guess.

In the frequency domain, one of the most reliable methodologies is the harmonic balance. In [START_REF] Lu | A modified incremental harmonic balance method for rotary periodic motions[END_REF], an extension of the harmonic balance method for computing the rotary and oscillatory periodic motion of a nonlinear smooth system is proposed. Mixed timefrequency-domain approaches are presented in [START_REF] Bonnin | A mixed timefrequency-domain approach for the analysis of a hysteretic oscillator[END_REF][START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF]. In the first paper, a nonlinear oscillator is analysed by linearizing the system along the solution predicted by the harmonic balance technique and then by computing the Floquet's multipliers by using a time-domain numerical algorithm. Instead, in [START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF], the harmonic balance method is implemented together with the envelope following method in time domain. Such approach is used to compute the steady-state behaviour and the associated period of nonlinear circuits forced by two input signals with different oscillation frequencies. If one is interested in approximate results, the harmonic balance with a single harmonic leads to the describing function method [START_REF] Duarte | Fractional describing function of systems with Coulomb friction[END_REF][START_REF] Huang | Steady-state analysis for a class of sliding mode controlled systems using describing function method[END_REF]. Such technique provides simple results about the existence of periodic solution and its parameters (the amplitude and the period), but it may fail, particularly when the system under consideration does not filter out the higher-order harmonics [START_REF] Engelberg | Limitations of the describing function for limit cycle prediction[END_REF].

Complementarity models have shown to be useful for representing PWL systems [START_REF] Vasca | A new perspective for modeling power electronics converters: complementarity framework[END_REF][START_REF] Schumacher | Complementarity systems in optimization[END_REF][START_REF] Sessa | A complementarity model for closed-loop power converters[END_REF][START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF][START_REF] Heemels | Equivalence of hybrid dynamical model[END_REF] and for computing their steady-state solutions [START_REF] Iannelli | Computation of limit cycles and forced oscillations in discrete-time piecewise linear feedback systems through a complementarity approach[END_REF][START_REF] Iannelli | Computation of limit cycles in Lur'e systems[END_REF][START_REF] Sessa | Mixed linear complementarity problems for the analysis of limit cycles in piecewise linear systems[END_REF]. In [START_REF] Iannelli | Computation of limit cycles and forced oscillations in discrete-time piecewise linear feedback systems through a complementarity approach[END_REF], the authors proposed the use of a linear complementarity approach for the evaluation of permanent oscillations in relay feedback systems with periodic exogenuous signals. The autonomous case for the particular class of Lur'e systems has been analysed in [START_REF] Iannelli | Computation of limit cycles in Lur'e systems[END_REF] by representing the discretized closed-loop system in a linear complementarity form. In [START_REF] Sessa | Mixed linear complementarity problems for the analysis of limit cycles in piecewise linear systems[END_REF], the periodicity condition was included as a constraint of the complementarity problem, but an a priori estimation of the period was still required. Instead, in [START_REF] Sessa | Computing period and shape of oscillations in piecewise linear Lur'e systems: a complementarity approach[END_REF], the period is considered as a further unknown and it is computed together with the periodic solution by constructing a suitable mixed complementarity problem. However, such approach fails if the solution does not satisfy the assumption of being sufficiently smooth.

In this paper, we propose the use of mixed quadratic complementarity problems (MQCPs) for the computation of periodic solutions exhibited by PWL systems. The analysis here proposed extends the approach presented in [START_REF] Sessa | Computing period and shape of oscillations in piecewise linear Lur'e systems: a complementarity approach[END_REF] by allowing nonsmooth solutions and by considering a more general class of PWL systems which includes those analysed in [START_REF] Iannelli | Computation of limit cycles and forced oscillations in discrete-time piecewise linear feedback systems through a complementarity approach[END_REF][START_REF] Iannelli | Computation of limit cycles in Lur'e systems[END_REF][START_REF] Sessa | Mixed linear complementarity problems for the analysis of limit cycles in piecewise linear systems[END_REF][START_REF] Sessa | Computing period and shape of oscillations in piecewise linear Lur'e systems: a complementarity approach[END_REF]. In particular, the PWL system is characterized by a nonlinearity which is representable as a linear combination of saturations and relays characteristics. The effectiveness of the proposed technique is shown by investigating several examples whose steadystate solutions cannot be obtained by using the complementarity approaches previously presented in the literature.

The paper is organized as follows. In Sect. 2, we present the autonomous PWL systems of interest and we give some preliminaries about complementarity problems and the complementarity representation of a class of PWL characteristics. Section 3 shows how to formulate a mixed quadratic complementarity problem (MQCP) for computing a periodic solution together with the period. In Sect. 4, two specific cases are analysed: when the PWL system has multiple (constant and/or nonconstant) periodic solutions, and when the period of the oscillation is known. In Sect. 5, we show the effectiveness of the proposed approach by considering different steady-state behaviours in several PWL systems of practical interest. In particular, we consider a stable periodic solution in a neural oscillator, an unstable periodic solution in a deadzone feedback system, a sliding periodic solution in a stick-slip system, the periodic oscillation belonging to a repressilator and a grazing-sliding bifurcation in a relay feedback system. The paper is concluded in Sect. 6.

Preliminaries

Let us introduce first some notation and preliminary definitions. Let vec(•) indicate the vectorization operator, 1 R N + is the set of nonnegative N -dimensional real vectors, I N denotes the N × N identity matrix and 1 N , ∞ N are the N -dimensional vectors whose components are ones and infinity, respectively. The symbol ⊗ indicates the Kronecker product. The following block circulant matrix notation is adopted

N (B 1 , B 2 ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ B 1 0 • • • • • • 0 B 2 B 2 B 1 0 • • • • • • 0 0 B 2 B 1 0 • • • . . . . . . . . . . . . . . . . . . 0 • • • 0 B 2 B 1 0 0 0 • • • 0 B 2 B 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 1
)
where B 1 and B 2 are square matrices of the same dimension and N is the number of times that B 1 is repeated on the main diagonal.

The mixed quadratic complementarity problem (MQCP) is another important ingredient of this paper. It is defined as follows. Given a function ϕ(z) : R r → R r , a square matrix M ∈ R r ×r , a vector q ∈ R r , and lower and upper bounds ,

u ∈ R r ∪ {-∞, +∞} r , a MQCP is to find z ∈ R r , w ∈ R r + , v ∈ R r + (2a) s.t. ϕ(z) + Mz + q = w -v (2b) z u, (z -) w = 0, (u -z) v = 0, ( 2c
)
where ϕ(z) is a vector of quadratic forms in z and the inequalities are meant componentwise. If and u are finite, the continuity of F(z) = ϕ(z) + Mz + q ensures the existence of solutions. More results about the existence of solutions for complementarity problems and variational inequality (VI) can be found in [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. Recall 1 vec{A} denotes the column vector obtained by stacking in one column all the columns of the matrix A. By vec(c 1 , c 2 , . . . , c n ) we will denote the column obtained by stacking all the c i columns, even when they have a different number of components.

that every (nonlinear) complementarity problem as in ( 2) is a variational inequality (VI). When the term ϕ(z) does not appear in (2a), the MQCP becomes a mixed linear complementarity problem (MLCP) [START_REF] Dirkse | The PATH solver: a nonmonotone stabilization scheme for mixed complementarity problems[END_REF]. Moreover, if the upper bound u has all infinity components and = 0, one gets v = 0 and the problem reduces to the classical linear complementarity problem with w and z being the usual complementarity variables [START_REF] Cottle | The Linear Complementarity Problem[END_REF].

The class of systems of interest can be represented by the block scheme in Fig. 1, where the dynamical system is described by

ẋ = A u x + B u u + e u (3a) ξ = C u x + D u u + f u (3b) u = -H u λ (3c) y = H y ξ, ( 3d
) with A u ∈ R n×n , B u ∈ R n× p , C u ∈ R p×n , D u ∈ R p× p , e u ∈ R n , f u ∈ R p , y = vec(y 1 , y 2 , . . . , y m ), λ = vec(λ 1 , λ 2 , . . . , λ m ), H y ∈ R m× p and H u ∈ R p×m .
Each block R i consists of a scalar saturationlike characteristic as in Fig. 2. Such characteristic with finite β R i , i , u i ∈ R, i < u i and μ i 0 can be represented in the following mixed linear complementarity form

μ i λ i -y i + β R i = w -v (4a) i λ i u i , (λ i -i )w = 0, (u i -λ i )v = 0, (4b)
with w 0 and v 0. If μ i > 0, for any given y i ∈ R the MLCP (4) has a unique solution λ i [START_REF] Gowda | Stability analysis of variational inequalities and nonlinear complementarity problems, via the mixed linear complementarity problem and degree theory[END_REF]. Therefore, in order to show that the MLCP (4) corresponds to the relation in Fig. 2, it is enough to verify that the feasible values of λ i correspond to values of y i belonging to the saturation characteristic. In particular, if λ i = i from (4b) it follows v = 0 and, since w 0, one obtains y i μ i i + β R i , that corresponds to the lower constant saturated piece of the characteristic. If i < λ i < u i from the second and third constraints in (4b) it follows w = v = 0, and from (4a) one obtains y i = μ i λ i + β R i that corresponds to the linear interval of the saturation characteristic. Finally, if λ i = u i , one obtains w = 0 and y i μ i u i +β R i , which corresponds to the upper saturation of the characteristic. If μ i = 0 the model (4) represents set-valued characteristics. Indeed, it is easy to verify that for μ i = 0 the model (4) represents a set-valued step function.

y 1 ym λ 1 λ 2 . . . λm - y i λ i i • 1 µ i • β R i μ i i + β R i μ i u i u i + β R i
Fig. 2 PWL saturation-like characteristic with μ i > 0

By collecting (4) and by defining R = vec(R 1 , R 2 , . . . , R m ), one can write the relation between y and λ as follows

find λ ∈ R m , w ∈ R m + , v ∈ R m + (5a) s.t. M R λ -y + β R = w -v (5b) R λ u R , (λ -R ) w = 0, (u R -λ) v = 0, (5c)
with

M R = diag{μ 1 , μ 2 , . . . , μ r } ∈ R m×m , R = vec( 1 , 2 , . . . , m ) ∈ R m , u = vec(u 1 , u 2 , . . . , u m ) ∈ R m and β R = vec(β R 1 , β R 2 , . . . , β R m ) ∈ R m .
The model ( 5) includes the representation of many typical PWL functions such as minimum, maximum, relay, deadzone and PWL characteristics obtained by linear combinations of these ones.

We now obtain a compact representation for the closed-loop system in Fig. 2. By substituting (3c) in (3a) and (3b), and (3b) in (3d) and by taking

A = A u (6a) B = B u H u (6b) C = H y C u (6c) D = H y D u H u (6d) e = e u (6e) f = H y f u , ( 6f
)
we obtain the following dynamical system ẋ = Ax + B(-λ) + e (7a)

y = C x + D(-λ) + f (7b) (y, λ) ∈ R, ( 7c
)
with R = vec(R 1 , R 2 , . . . , R m ) given by ( 5) and

(y i , λ i ) ∈ R i , i = 1, . . . , m.
The class of systems [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] is the one considered in this paper where (A, B, C) is a minimal state space real-

ization with A ∈ R n×n , B ∈ R n×m , C ∈ R m×n , D ∈ R m×m , e ∈ R n
and f ∈ R m all being constant and the time derivative is meant almost everywhere. It is assumed that for every initial condition x(t 0 ) the system [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] has an absolutely continuous solution x(t) : [t 0 , +∞) → R n that satisfies [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] for almost every t t 0 . A solution x(t) of the system ( 7) is periodic if there exists a positive real number T such that x(t + T ) = x(t) for any t 0. The stability concepts and definitions typically used for systems and equilibria can be also applied to periodic solutions which can be stable, asymptotically stable and unstable [START_REF] Farkas | Periodic Motions. Series in Applied Mathematical[END_REF].

In the next section, it is shown how the combination of [START_REF] Melin | Two types of limit cycles of a resonant converter modelled by threedimensional system[END_REF] with the discretized version of (7a)-(7b) allows to construct a MQCP for the computation of periodic solutions of [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF].

Computation of periodic oscillations

Let x(t) be a nonconstant periodic solution of [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] with unknown period, say T . The dynamical model ( 7) can be normalized with respect to the unknown period by using the time scaling t = T τ , where τ is a dimensionless time variable. Then (7) can be rewritten as

x = T Ax + T B(-λ) + T e (8a) y = C x + D(-λ) + f (8b) (y, λ) ∈ R, ( 8c
)
where x is the derivative with respect to τ . Any periodic solution of (7) with period T will correspond to a periodic solution of ( 8) with period 1. Note that the right-hand side of (8a) is quadratic with respect to the unknowns T, x and λ. By using ( 5) and ( 8), the problem of computing the periodic solution is to

find λ(τ ) : [0, 1] → R m , w(τ ) : [0, 1] → R m + , v(τ) : [0, 1] → R m + , T ∈ R, x(τ ) : [0, 1] → R n (9a) s. t. T Bλ(τ ) -T Ax(τ ) + x (τ ) -T e = 0 (9b) (M R + D)λ(τ ) -C x(τ ) -f + β R = w(τ ) -v(τ ) (9c) R λ(τ ) u R , (λ(τ) -R ) w(τ ) = 0, (u R -λ(τ )) v(τ ) = 0 ( 9 d ) λ(0) = λ(1), x(0) = x(1), w(0) = w(1), v(0) = v(1), ( 9e
)
where the constraints (9b)-(9d) must hold for any τ ∈

[0, 1].
The phase of a periodic solution of an autonomous system is not fixed. Then any time translation of a periodic solution provides another 'different' periodic solution. In other words, if the PWL system admits a periodic solution, it admits an infinite number of periodic solutions each one differing from the others by a translation in time. In order to fix the initial phase of the periodic solution, one more equation is required, which is usually called anchor equation [START_REF] Doedel | AUTO: a program for the automatic bifurcation analysis of autonomous systems[END_REF][START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. A possible anchor equation is the phase condition proposed in [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]:

x j ( τ ) = T c j Ax( τ ) -Bλ( τ ) + e = 0, (10)
where x j is a generic j-th component of the state at an arbitrarily chosen time instant τ ∈ [0, 1] and c j ∈ R n is a vector with all elements equal to zero except for the j-th element equal to 1. Note that the index j can be chosen arbitrarily for sufficiently smooth solutions. Indeed, in case of periodic solutions x(t) of class C 1 , the time derivative of each state variable must be zero at least at one time instant τ ∈ [0, 1].

In order to solve ( 9)- [START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF], one can discretize (9b) and reformulate the problem in terms of a discretetime complementarity problem. By using the (θ, γ ) discretization technique [START_REF] Acary | Time-stepping via complementarity[END_REF] with a sampling period 1/N , N being an integer, and by using the subscript k for indicating the k-th sample of a variable, from (9b) one obtains

T γ Bλ k +T (1-γ )Bλ k-1 -T θ Ax k -T (1-θ)Ax k-1 + N x k -N x k-1 -T e = 0, (11)
where θ ∈ [0, 1] and γ ∈ [0, 1]. The constraints (9c)-(9d) can be written for each sampling time instant:

(M R + D)λ k -C x k -f + β R = w k -v k (12a) R λ k u R , (λ k -R ) w k = 0, (u R -λ k ) v k = 0, (12b)
with w k 0, v k 0 and k = 1, . . . , N . Moreover, (9e) can be rewritten as

λ 0 = λ N , x 0 = x N , w 0 = w N , v 0 = v N . ( 13
)
Finally, one can assume that τ in ( 10) is a sampling time instant. Therefore, by defining k = τ N , the anchor equation ( 10) can be written as

T c j Ax k -Bλ k + e = 0, (14)
or, equivalently,

x k -x k-1 = 0. ( 15
)
By collecting ( 11)-( 13) for k = 1, . . . , N in a matrix form we are now able to formulate the problem of computing the periodic solution (and its period) as a MQCP. In particular, the unknown vector z for the MQCP is given by the samples λ k and x k for k = 1, . . . , N and the period T :

z = vec(λ, x, T ) = vec(λ 1 , . . . , λ N , x 1 , . . . , x N , T ). ( 16
)
The lower and upper bounds for T and for each x k with k = 1, . . . , N are set as -∞ and +∞:

= vec( λ , x , l T ) = vec(1 N ⊗ R , -∞ N •n , -∞) (17a) u = vec(u λ , u x , u T ) = vec(1 N ⊗ u R , +∞ N •n , +∞). ( 17b
)
Given the discretization parameters (θ, γ ) and the number of discrete samples per period N , by considering ( 16)-( 17), by collecting ( 11)-( 12) for k = 1, . . . , N together with [START_REF] Leine | Stick-slip vibrations induced by alternate friction models[END_REF] and by considering ( 14), the MQCP can be represented in the form [START_REF] Camlibel | Switched networks and complementarity[END_REF] with

M = ⎡ ⎣ 0 G e M R C 0 h λ h x 0 ⎤ ⎦ ∈ R (N n+N m+1)×(N n+N m+1) (18a) q = vec(0 N n , f , δ q ) ∈ R N n+N m+1 (18b) ϕ(z) = vec(φ 1 (z), . . . , φ N (z), 0 N m+1 ), ( 18c
)
where φ k (z) : R N n+N m+1 → R n are given by

φ k (z) = T γ Bλ k + T (1 -γ )Bλ k-1 -T θ Ax k -T (1 -θ)Ax k-1 , (19)
for k = 1, . . . , N together with [START_REF] Leine | Stick-slip vibrations induced by alternate friction models[END_REF], and

G = N N (I n , -I n ) (20a) e = -1 N ⊗ e (20b) M R = I N ⊗ (M R + D) (20c) C = -I N ⊗ C (20d) f = 1 N ⊗ (-f + β R ), ( 20e
)
h λ ∈ R N m is a vector with all zero elements, h x ∈ R N n is a vector with all zero entries but two elements -1 and +1 corresponding to the j-th component of x k-1 and x k , representing the anchor equation ( 15), and δ q = 0. The described formulation is used to compute a stable periodic oscillation in a neural oscillator and an unstable periodic solution in a deadzone feedback system in Sects. 5.1 and 5.2, respectively. Indeed, in such examples the feedback characteristic is Lipschitz continuous, then a solution of class C 1 is expected. If x(t) is of class C 0 , as in the case of a relay feedback system, we have to pay attention to the fact that the time derivative can jump, so a different phase constraint should be derived. In particular, we assume that there exists a component of the system output, say c j y, which is equal to a value α j at a time instant τ ∈ [0, 1]:

c j y = c j C x( τ ) -Dλ( τ ) + f = α j . ( 21
)
The constant α j should be selected by exploiting the system structure. In particular, since we use the phase condition (21) because the trajectory x(t) is only C 0 , the discontinuity of the time derivative ẋ(t) allows to say that there exist some j and τ ∈ [0, 1] such that (21) holds with α j being a discontinuity point of the relation R j .

As in the previous case, we write [START_REF] Engelberg | Limitations of the describing function for limit cycle prediction[END_REF] in discrete time as follows

c j C x k -Dλ k + f = α j , ( 22
)
with k = τ N . The anchor equation ( 22) can be included in the form (18) by setting h λ and h x with all zero entries except for the components corresponding to k which are equal to -c j D and c j C, respectively, and δ q = c j f -α j . In Sects. 5.3, 5.4 and 5.5, the anchor equation ( 22) is used to compute periodic nonsmooth oscillations in a stick-slip system, in a repressilator and in a relay feedback system.

The continuous PWL interpolation of the sequence x k obtained by solving the MQCP derived above is expected to approximate better and better the continuous-time periodic solution x(t) of ( 7) as the number of samples per period N increases. Such convergence property of the discrete-time solution is also called 'consistency' [START_REF] Han | Convergence of time-stepping for passive and extended linear complementarity systems[END_REF]. A formal proof of consistency of the solution of time-stepping methods used for solving boundary-value differential variational inequalities can be found in [START_REF] Pang | Differential variational inequalities[END_REF]. However, the nonsingularity condition of matrices involved in the boundary-value constraint is not valid in the case of periodic solutions.

Multiple solutions and known period

In this section, the MQCP formulated above is specified under two interesting conditions, i.e. when more periodic solutions coexist and when the period of the oscillation is known.

Elimination of a given periodic solution

Depending on the relation (y, λ) ∈ R, the corresponding PWL system can present multiple periodic solutions as well as constant solutions. For instance, when e and f are zero and (0, 0) ∈ R, the origin is also a solution (an equilibrium) of the system. Once a periodic (constant or nonconstant) solution is computed, say x, one may be interested in finding another periodic trajectory. To this aim, one can formulate an analogous MQCP by adding a further constraint which excludes the solution x from its feasibility set. In particular, the following condition can be used:

x (•, î) -x (•, î) x (•, î) -x (•, î) > 0 ( 23
)
for some chosen î, where x (•, î) ∈ R N is the vector obtained by collecting all samples of the î-th component of the state vector. The condition [START_REF] Schumacher | Complementarity systems in optimization[END_REF] can be included in the MQCP by using the following complementarity representation:

find z ρ ∈ R, w ρ ∈ R + , v ρ ∈ R + (24a) s. t. x (•, î) -x (•, î) x (•, î) -x (•, î) -= w ρ -v ρ (24b) ρ z ρ u ρ , (z ρ -ρ )w ρ = 0, (u ρ -z ρ )v ρ = 0, (24c)
where is a small positive parameter, with ρ = 0 and u ρ equal to +∞ (so that v ρ = 0). Indeed, from the nonnegativity of the variable w ρ , the constraints (24b)

with v ρ = 0 implies that (x (•, î) -x (•, î) ) (x (•, î) - x (•, î)
) must be strictly positive; i.e. there exists at least a sample of x (•, î) which is different from a sample of the computed current solution x (•, î) .

MLCP for known period

A particular case of MQCP occurs when the period T of the periodic solution is known. In this case, the periodic solution computation can be formulated in terms of a MLCP. Since T is known, Eqs. ( 11) are linear with respect to the unknowns λ k and x k . By collecting [START_REF] Theodosiou | On periodic steady state response and stability of Filippov-type mechanical models[END_REF] and (12a) for k = 1, . . . , N , by defining and u as in [START_REF] Bonnin | A mixed timefrequency-domain approach for the analysis of a hysteretic oscillator[END_REF] without the last element and by using the periodicity condition [START_REF] Leine | Stick-slip vibrations induced by alternate friction models[END_REF], one obtains

B λ + Ax + e = 0 (25a) M R λ + C x + f = w -v, ( 25b
)
with

A = N I n - T N θ A, -I n - T N (1 -θ)A (26a) B = T N N (γ B, (1 -γ )B) , (26b)
and e, M R , C and f defined as in (20b)-(20e), respectively.

Since T is known, the anchor equation ( 10) is not required anymore. In this case, we can define z = col(λ, x) and formulate the above problem as a MLCP, that is (2) without ϕ(z), with the following matrices

M = M R C B A and q = f e . ( 27
)
In the particular case of a nonsingular A, one can solve (25a) for the vector x thus obtaining

x = -A -1 (B λ + e). ( 28
)
By substituting [START_REF] Iannelli | Computation of limit cycles in Lur'e systems[END_REF] in (25b) one obtains

M R -C A -1 B λ + f -C A -1 e = w -v, ( 29
)
which is in the form (2b) without ϕ(z) and with z = λ, M = M R -C A -1 B and q = f -C A -1 e. The solution λ of the corresponding MLCP can be used for the computation of the periodic solution x by using [START_REF] Iannelli | Computation of limit cycles in Lur'e systems[END_REF].

Note that for T /N → 0, the determinant of A tends to zero. Then we can state that when a high accuracy is required, i.e. a small step size T /N , the formulation in ( 25) is preferable to [START_REF] Sessa | Mixed linear complementarity problems for the analysis of limit cycles in piecewise linear systems[END_REF].

In the case that the PWL system presents constant solutions, the complementarity problem can be reformulated so that the new problem has no constant solutions. Since the use of ( 24) makes the problem nonlinear, in order to preserve the linearity we can use an alternative approach.

As an example, let us assume that the PWL system has the origin as an equilibrium point, as discussed at the beginning of Sect. 4.1. In order to find only nonzero solutions, one can impose (without loss of generality) x k to be greater than some positive value ε for some chosen k. In particular, the constraint x k ε can be written in the complementarity form by modifying the constraints on some samples of the vector x. For instance, one can set x k,1 = ε and u x k,1 = +∞, where x k,1 , u x k,1 are the lower bound and the upper bound of the first component of x at the time k. Note that the approach presented above is not required in the case of an unknown period; indeed, such further constraint affects the phase condition of the solution which is already determined by the anchor equation.

Examples

In this section, the effectiveness of the proposed technique for computing period and waveform of oscillations is shown by considering PWL systems with different steady-state behaviours. We present a neural oscillator and a deadzone feedback system, which have solutions belonging to class C 1 . The second group of examples, including a stick-slip system, a repressilator and a relay feedback system, present oscillations of class C 0 . The computation of the periodic solution and the corresponding period is obtained by solving a MQCP with the matrices given in [START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF].

For all the following examples, we have computed the 'exact' steady-state solution of the discretized system by constructing the nonlinear closed-loop map. The idea for obtaining such map is to consider the solutions of the differential state equations, which model the 'modes' of the system, and to cascade them in order to link a sampled state to the next sample. It is important to highlight that the exact solution can be analytically obtained only if the sequence of modes is known, while the proposed approach is formulated without considering any hypothesis of such sequence. We numerically show that by varying the number of samples per period, i.e. by improving the resolution of the discretization, the error between the solution obtained with the MQCP and the exact one decreases. The maximum number of samples is chosen for a reasonable computational effort to perform the MATLAB code on a Intel Core i7 clocked at 2.40 GHz. The complementarity problems have been solved by using the PATH solver [START_REF] Dirkse | The PATH solver: a nonmonotone stabilization scheme for mixed complementarity problems[END_REF], that is a particular version of the Newton method for nonsmooth problems. The solver requires an initial guess, which by default is set to zero. Since PATH provides only a local convergence, it is clear that if the initial point is close to an existing solution, one has more chances that it converges to that particular solution; e.g. one can choose as initial guess for the period the one predicted by the describing function. In the following examples, we report the initial point used in the PATH solver, whose value is chosen from a set of trivial vectors and is not related to neither to the searched solution nor to the corresponding sequence of modes. The numerical results are also compared by varying the discretization parameters.

Stable periodic solution in a neural oscillator

The dynamical model of a neural network consisting of η mutually inhibiting neurons with adaptation can be represented in the following form [START_REF] Matsuoka | Sustained oscillations generated by mutually inhibiting neurons with adaptation[END_REF]:

τ 1 χi = -χ i -κv i -ν η j=1, j =i u j + e c ( 30a
)
τ 2 vi = -v i + u i (30b) u i = max{0, χ i }, (30c)
with i = 1, . . . , η. Here, χ i is the membrane potential of the i-th neuron, v i represents the degree of adaptation, e c is the total input from outside the network that is assumed positive and constant with the time, τ 1 is a time constant, τ 2 and κ are the parameters that specify the time course of the adaptation, and ν indicates the strength of the inhibitory connection between the neurons. In [START_REF] Matsuoka | Sustained oscillations generated by mutually inhibiting neurons with adaptation[END_REF], it is shown that for certain values of model parameters, the neural network has no stable equilibrium points and produces sustained oscillations. In the following, we consider a neural network with two neurons; i.e. η = 2, see Fig. 3. For neuron 1 and neuron 2, let us denote by x 1 and x 3 the state variables that represent the membrane potential, respectively, and by x 2 and x 4 the state variables for the degree of selfinhibition, respectively. Then the neural oscillator can be represented in the form (3) with the following matrices:

A u = ⎡ ⎢ ⎢ ⎢ ⎣ -1 τ 1 -κ τ 1 0 0 0 -1 τ 2 0 0 0 0 -1 τ 1 -κ τ 1 0 0 0 -1 τ 2 ⎤ ⎥ ⎥ ⎥ ⎦ , B u = ⎡ ⎢ ⎢ ⎢ ⎣ 0 -ν τ 1 1 τ 2 0 -ν τ 1 0 0 1 τ 2 ⎤ ⎥ ⎥ ⎥ ⎦ , e u = ⎡ ⎢ ⎢ ⎣ e c 0 e c 0 ⎤ ⎥ ⎥ ⎦ , ( 31a
) C u = 1 0 0 0 0 0 1 0 , D u = 0 0 0 0 , f u = 0 0 . ( 31b
)
The relation R = vec(R 1 , R 2 ) represents the inhibitory connection between the two neurons, with R i that corresponds to λ i = max{0, y i }, i = 1, 2, so as indicated by (30c). Each R i can be modelled as in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] with μ i = 1, β R i = 0, i = 0 and u i = +∞, and they can be 'linked' to the system by choosing H u = -I 2 and H y = I 2 , in (3c) and (3d), respectively. By using ( 6), the neural model network is recast in form (7a) and (7b) and the feedback relation is In [START_REF] Goncalves | Region of stability for limit cycles in piecewise linear systems[END_REF], it is shown that by considering the parameters τ 1 = 0.1 s, τ 2 = 0.2 s, κ = 2, ν = 2 and e c = 1, the system has a locally stable periodic oscillation. Figure 4 shows the solution computed by solving the corresponding MQCP with N = 600, θ = 0.5 and γ = 0.5 that yields a value of the period T = 0.8973 s. The neural system has also a constant unstable solution in [0.2, 0.2, 0.2, 0.2] that was eliminated by following the procedure presented in Sect. 4.1. The initial guess for the solver has been set as z 0 = [x 0 ; T 0 ; λ 0 ; z 0 ρ ], where x 0 = [1, . . . , 1, -1, . . . , -1], λ 0 = [1, . . . , 1, 0, . . . , 0, 1, . . . , 1], T 0 = 0.897 s, and z 0 ρ = 0. Let us consider the error between the exact solution and the solution computed with the complementarity approach by varying the number of samples per period and the parameters θ and γ . Since the relation R is a Lipschitz continuous function of x, the solution is expected to belong to the class C 1 . As suggested in [START_REF] Acary | Time-stepping numerical simulation of switched circuits within the nonsmooth dynamical systems approach[END_REF], in the case of solutions of class C 1 , it is possible to compute the solution of the discretized system by varying both θ and γ in order to obtain the pair which provides the best accuracy for the solution. Figure 5 shows the evolution of the maximum error in the computation of the x 2 trajectory, while in Fig. 6 it reported the error for the period computation.

As one could expect, the error decreases when the number of samples increases. For θ = 0.5 and γ assuming values between 0.5 and 1, the accuracy is improved, but we do not achieve order 2 when γ = 0.5. This in Fig. 6 Neural system. Error in the period computation obtained with different discretization techniques and with N ∈ [START_REF] Bernardo | Piecewise-Smooth Dynamical Systems[END_REF]600] mainly due to the lack of regularity of the computed solution. This phenomenon is generally observed when higher-order methods are used for a solution with limited smoothness; see [4, §9.1].

Unstable periodic solution in a deadzone feedback system

The deadzone behaviour is typical in many real actuator systems, including mechanical connections, hydraulic y λ 1 -1 Fig. 7 Deadzone characteristic servo valves, piezoelectric translators and electronic circuits. A graphical representation is given in Fig. 7, and an analytical description is presented below

λ = ⎧ ⎪ ⎨ ⎪ ⎩ y -1 if y > 1 0 i f -1 y 1 y + 1 if y < -1 . ( 32
)
In this section, we analyse the steady-state periodic behaviour of the dynamical system in form (3) with the following matrices:

A u = ⎡ ⎣ 0 1 0 0 0 1 0 -3 -4 ⎤ ⎦ , B u = ⎡ ⎣ 0 0 20 ⎤ ⎦ , e u = ⎡ ⎣ 0 0 0 ⎤ ⎦ , (33a) C u = 1 0 0 , D u = 0, f u = 0, (33b)
and whose feedback characteristic is the deadzone in Fig. 7. Note that A is not Hurwitz, then self-induced oscillation is possible [START_REF] Hu | A unified Lyapunov approach to analysis of oscillations and stability for systems with piecewise linear elements[END_REF].

The dynamical system is a low-pass filter and the polar diagram intersects the negative real axis, so the closed-loop system is a good candidate for the application of the describing function technique [START_REF] Khalil | Nonlinear Systems[END_REF]. The value of the approximated period is computed by considering the oscillation frequency ω * corresponding to the point in which the polar diagram of G( jω) (continuous line in Fig. 8) intersects the negative reciprocal of the describing function graph (dashed line in Fig. 8). We obtain ω * = 1.723 rad/s, then the corresponding period T * = 3.627 s. Moreover, the orientation of the two curves predicts that the periodic oscillation is unstable. The following numerical results show that the proposed complementarity approach is able to get the periodic behaviour together with the period also in the case of unstable solutions. The deadzone characteristic can be expressed as the difference of two saturation-like characteristics, R 1 and R 2 . In particular, R 1 corresponds to λ 1 = y 1 , while R 2 is λ 2 = sat(y 2 ), that is the unitary saturation function. These two relations can be modelled as in (4) by choosing for R 1 , μ 1 = 0, β R 1 = 0, 1 = -∞ and u 1 = +∞, while for R 2 , μ 2 = 1, β R 2 = 0, 2 = -1 and u 2 = +1. Then we choose H u = 1, -1 and H y = vec(1, 1) in (3c) and (3d), respectively. By using [START_REF] Yang | Limit cycle and chaotic invariant sets in autonomous hybrid planar systems[END_REF], the dynamical system can be recast in form (7a) and (7b), together with the complementarity representation of feedback relation discussed above.

10

Figure 9 shows the periodic solution computed through the proposed technique, by solving a MQCP with matrices in [START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF] and, with N = 5400, θ = 0.5 and γ = 0.5, which yields to a value of the period T = 3.6620 s. The initial guess for the solver has been set as z 0 = [x 0 ; T 0 ; λ 0 ; z 0 ρ ], where x 0 = [1, . . . , 1, -1, . . . , -1], λ 0 = [1, . . . , 1, -1, . . . , -1], T 0 = 3.6 s, and z 0 ρ = 0. The dotted line graph corresponds to the numerical result obtained by implementing a time-stepping simulation of the initial value problem [START_REF] Acary | Time-stepping numerical simulation of switched circuits within the nonsmooth dynamical systems approach[END_REF][START_REF] Heemels | A time-stepping method for relay systems[END_REF], with an initial condition sufficiently close to the periodic solution (x 0 = [1.3, -7.848, 4]). Of course the standard time-stepping simulation is not able to show the unstable orbit: the proposed approach is particularly advantageous in such cases.

The error between the exact solution and the one computed with the complementarity approach has been analysed by varying the number of samples per period between 100 and 5400 and by considering values of the parameters θ and γ in the interval [0.5, 1]; see Figs. 10 and11. It is possible to prove that the trivial solution is also a solution of the system considered in this example, so the constraint [START_REF] Sessa | A complementarity model for closed-loop power converters[END_REF] was added to the problem [START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF] to directly get the periodic solution. Qualitatively, the same comments presented for the previous example about the accuracy of the error can be repeated; see Figs. 5 and 6.

Sliding periodic solution in a stick-slip system

Let us consider a mass m connected to a spring with elastic constant K , which is pulled at a constant speed e s . Let x 1 be the elongation of the spring, and x 2 the mass velocity. Such system can be represented in the form [START_REF] Acary | Time-stepping numerical simulation of switched circuits within the nonsmooth dynamical systems approach[END_REF] with

A u = 0 -1 K m 0 , B u = 0 1 m , e u = e s 0 , (34a)
C u = 0 1 , D u = 0, f u = 0, (34b)
and R being the friction characteristic represented in Fig. 12, where F c and F s are the Coulomb friction and the stiction force, respectively. As in the previous example, such characteristic can be obtained as the combination of two saturation-like characteristics, say R 1 and R 2 . In this case, R 1 represents λ 1 = relay(y 1 ), where relay(•) is the unitary relay characteristic, and

R 2 is λ 2 = sat m R (y 2 )
, where sat m R (•) is the unitary saturation function where the nonzero slope of the characteristic is 1/m R . The two relations can be modelled as in (4) by choosing for

R 1 , μ 1 = 0, β R 1 = 0, 1 = -1 and u 1 = +1, and for R 2 , μ 2 = m R , β R 2 = 0, 2 = -1 and u 2 = +1. y λ -F c 1/m R -F s F s F c
Fig. 12 Coulomb friction characteristic with Stribeck effect with null viscous friction: y is the mass velocity and λ is the friction force

Then we choose

H u = [F s , -(F s -F c )] and H y = vec (1, 1) in (3c) and (3d)
, respectively, and we write the dynamical system in form [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] with the matrices given by [START_REF] Yang | Limit cycle and chaotic invariant sets in autonomous hybrid planar systems[END_REF].

Let us select m = 1 kg, K = 2 N/m, F c = 2.94 N, F s = 5.88 N, m R = 1/12, and e s = 2 m/s. In [START_REF] Astrom | Revisiting the LuGre friction model[END_REF], it is shown that this system can exhibit complex behaviours and that, depending on the selected integration technique, it is possible that sticking phases of the mass can be missed by the numerical simulation. Due to the discontinuity of the PWL relation at y = 0, a solution belonging to C 0 is expected. For such reason, the computation of the periodic solution and the corresponding period was obtained by solving the MQCP with matrices given in [START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF], using the anchor equation with the one in [START_REF] Vasca | A new perspective for modeling power electronics converters: complementarity framework[END_REF]. Figure 13 shows that the complementarity approach is able to get the steady-state periodic behaviour of the system, in which it is possible to recognize the stick and the slip modes. When the trajectory is in the stick mode, the friction force increases rapidly and stops the motion. The system remains in the stick mode as long as the spring force is smaller than the stiction force. After that, the system enters in the slip mode and the mass accelerates subject to the Coulomb friction. The spring is then compressed and the motion of the mass stops again and the spring force needs to win again the stiction force. Such solu- The initial guess for the solver has been set as z 0 = [x 0 ; T 0 ; λ 0 ; z 0 ρ ], where x 0 = [5, . . . , 5, 0, . . . , 0], λ 0 = [(F s -F c ), . . . , (F s -F c ), 0, . . . , 0], T 0 = 4.91 s, and z 0 ρ = 0. Since the smoothness of the solution is one degree lower with respect to the previous examples, according to the discretization scheme in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF], we have fixed γ = 1 and we have varied only θ. The evolution of the error shown in Figs. 14 and15 confirms the effectiveness of the proposed approach.

Nonsmooth stable periodic solution in a repressilator

Let us now consider the model of a repressilator, which is a synthetic genetic regulatory network that uses a cyclic repression structure. In the literature (see, e.g., [START_REF] Tonnelier | Cyclic negative feedback systems: what is the chance of oscillation?[END_REF] and the reference therein), it is shown the existence of oscillations for a three-dimensional system in which sigmoidal regulation functions that deter- mine the interaction between network elements are replaced by step functions. The resulting PWL systems have been proposed as a modelling framework in biology allowing efficient simulations, as confirmed in [START_REF] Acary | Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems[END_REF] where the complementarity approach is used to compute the behaviour of a repressilator in the bacterium Escherichia coli. Such a repressilator consists of three genes, and it can be modelled by considering the evolution of three variables which represent concentrations of the proteins Lacl, TetR and Cl. Such a repressilator model can be easily put in the form [START_REF] Acary | Time-stepping numerical simulation of switched circuits within the nonsmooth dynamical systems approach[END_REF] with

A u = ⎡ ⎣ -a 1 0 0 0 -a 2 0 0 0 -a 3 ⎤ ⎦ , B u = ⎡ ⎣ -b 1 0 0 0 -b 2 0 0 0 -b 3 ⎤ ⎦ , e u = ⎡ ⎣ e 1 e 2 e 3 ⎤ ⎦ , (35a) C u = ⎡ ⎣ 0 0 1 1 0 0 0 1 0 ⎤ ⎦ , D u = 0, f u = 0, (35b)
and

R = vec(R 1 , R 2 , R 3 )
where

λ i = ⎧ ⎨ ⎩ 0 ify i > 1 [0, 1] if y i = 1 1 ify i < 1 , ( 36
)
and i = 1, 2, 3. Each relation R i can be modelled as the MLCP (4) with μ i = 0, β R i = -1, i = 0 and u i = 1. Then we choose H u = I 3 and H y = -I 3 in (3c) and (3d), respectively, and we write the dynamical system in form [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] with the matrices given by ( 6). As showed in [START_REF] Acary | Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems[END_REF], the system has an unstable equilibrium point in (1, 1, 1) and a stable nonsmooth periodic solution. In order to compute such solution and the associated period, we solve a MQCP with matrices in [START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF] replacing the anchor equation with the one in [START_REF] Vasca | A new perspective for modeling power electronics converters: complementarity framework[END_REF]. The initial guess for the solver has been set as z 0 = [x 0 ; T 0 ; λ 0 ; z 0 ρ ], where x 0 = [3, . . . , 3, 0, . . . , 0], λ 0 = [0, 1, 0, 0, 1, 0, . . . , 0, 1, 0], T 0 = 14.1959 s, and z 0 ρ = 0. Figure 16 shows the time evolution of the periodic steady state for the state variables, obtained with the proposed approach. The error results in Figs. 17 and 18 show that, for this particular example, the solution is not much sensitive to the discretization parameters.

Grazing-sliding bifurcation and multiple periodic solutions in relay feedback systems

A relay feedback system may present a great variety of behaviours at steady state. In this section, we analyse two examples. In the first one, the system undergoes a grazing-sliding bifurcation, while the second example presents multiple periodic solutions; in particular, one solution is stable and another is unstable.

Let us consider the third-order system described by the following matrices The relay characteristic can be modelled as the MLCP (4) with M R = 0, β R = 0, = -1 and u = 1.

A u = ⎡ ⎣ -(2ζ ω + α) 1 0 -(2ζ ωα + ω 2 ) 0 1 -αω 2 0 0 ⎤ ⎦ , B u = ⎡ ⎣ κ 2κσρ κρ 2 ⎤ ⎦ , e u = ⎡ ⎣ 0 0 0 ⎤ ⎦ , ( 37a
) C u = 1 0 0 , D u = 0, f u = 0. ( 37b
We choose H u = 1 and H y = 1 in (3c) and (3d), respectively, and we write the dynamical system in form [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] with the matrices given by [START_REF] Yang | Limit cycle and chaotic invariant sets in autonomous hybrid planar systems[END_REF]. In [START_REF] Bernardo | Piecewise-Smooth Dynamical Systems[END_REF], it is shown that such system undergoes different type of bifurcations by varying the system parameters. We consider the grazing-sliding bifurcation that occurs when the parameter ζ varies between 0.082 and 0.098 with all the other parameters fixed at the values ω = 10.84, -σ = α = ρ = κ = 1. In Fig. 19, we report the periodic steady-state evolutions of the system output x 1 and of the corresponding relay output λ obtained with three different values of the parameter ζ . The grazing-sliding bifurcation is clearly visible. For ζ = 0.082, there is a stable, symmetric trajectory with 6 pieces of sliding, which is shown in Fig. 19a. The evolution of the relay output is presented in Fig. 19b and, as expected, the value of λ is between -1 and 1 during the six time intervals corresponding to the six sliding pieces. The trajectory undergoes a grazing-sliding bifurcation and becomes a trajectory with 4 pieces of sliding for ζ = 0.098 as shown in Fig. 19e,f. Figure 19c shows that for ζ = 0.091 the trajectory grazes the switching surface at time instants 0.3475 and 3.243 s. Indeed, at the same time instants, the corresponding value of λ is between -1 and 1, as shown in Fig. 19d. The initial point for the solver PATH is set as z 0 = [x 0 ; T 0 ; λ 0 ; z 0 ρ ], where x 0 = [1, . . . , 1, -1, . . . , -1], λ 0 = [-1, . . . , -1, 1, . . . , 1], T 0 = 8.642 s, and z 0 ρ = 0. Let us consider the relay feedback system with the following matrices

A u = ⎡ ⎣ -2 -1 0 1 0 0 0 1 0 ⎤ ⎦ , B u = ⎡ ⎣ 1 0 0 ⎤ ⎦ , e u = ⎡ ⎣ 0 0 0 ⎤ ⎦ , ( 38a
) C u = 0.02 0.3 1 , D u = 0, f u = 0. (38b)
As in Sect. 5.2, let us use the describing function technique for predicting the existence of periodic solutions. Figure 20 shows the Nyquist diagram of system [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] together with the negative reciprocal of the relay describing function, which corresponds to the negative real axis. The Nyquist diagram intersects the negative real axis in two points in which ω assumes the values 1.65 and 4 rad/s. This means that the closed-loop system has two periodic solutions with different period. In particular, the system has an unstable periodic solution of period T = 1.570 s and a stable periodic solution of period T = 3.808 s. By using the complementarity representation for the relay characteristic given above, we can write the dynamical system in form [START_REF] Aprille | A computer algorithm to determine the steady-state response of nonlinear oscillators[END_REF] with the matrices given by [START_REF] Yang | Limit cycle and chaotic invariant sets in autonomous hybrid planar systems[END_REF].

The unstable solution in Fig. 21 and the corresponding period are computed by solving the MQCP in [START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF] and by setting N = 400, γ = 1 and θ = 1. The computed value of the period is 1.8030 s. The initial point for the solver PATH is set as z 0 = [x 0 ; T 0 ; λ 0 ], where x 0 = [0, . . . , 0, 1, . . . , 1], λ 0 = [1, . . . , 1, 0, . . . , 0], and T 0 = 2 s. For the dynamical system under consideration, the constant trivial solution is also a solution. By using the arguments in Sect. 4, suitable constraints were added to the MQCP in order to eliminate both the computed periodic unstable solution and the trivial one. Then the periodic stable steady-state solution in Fig. 22 was obtained, whose period is 3.5645 s. We use the same initial point, except that we added the initial value for the variables z ρ 1 and z ρ 2 , which was set to zero.

Conclusions

The mixed quadratic complementarity approach has been proposed for the computation of periodic solu- Fig. [START_REF] Engelberg | Limitations of the describing function for limit cycle prediction[END_REF] Unstable periodic steady-state solution of the relay feedback system in [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. The solution is obtained by fixing N = 400, θ = 1 and γ = 1. The corresponding value of the period is 1.8030 s tions in a class of piecewise linear (PWL) systems, which consists of linear time-invariant systems with PWL feedback relations representable as linear combinations of saturation-like or step characteristics. For many practical systems belonging to this class, the period and the shape of the oscillation are difficult to be predicted, then phase conditions acting as anchor equations for the periodic solution have been added to the complementarity problem. It has been showed how to build a mixed quadratic complementarity problem (MQCP) which models the discretized closed-loop system together with the periodicity constraint and the Fig. [START_REF] Vasca | A new perspective for modeling power electronics converters: complementarity framework[END_REF] Stable periodic steady-state solution of the relay feedback system in [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. The solution is obtained by fixing N = 400, θ = 1 and γ = 1. The corresponding value of the period is 3.5645 s phase condition. The solution of such MQCP gives, as outcome, the periodic solution and the corresponding period. The effectiveness of the proposed approach has been tested by considering several practical PWL systems with different steady-state periodic behaviours: a stable periodic solution in a neural oscillator, an unstable oscillation for a deadzone feedback system, a sliding orbit in a stick-slip system, a periodic solution belonging to a repressilator, a grazing-sliding bifurcation in a relay feedback system and coexisting periodic solutions in a relay feedback system. The simulations, implemented by using different discretization schemes and sampling frequencies, showed the good accuracy of the proposed approach with respect to the analytical solution, computed by assuming the a priori knowledge of the modes sequence. Future work will investigate the use of the proposed tool for the computation of periodic solutions in a continuation framework.

Fig. 1

 1 Fig. 1 Block diagram of a PWL system ẋ = Aux + Buu + eu ξ = Cux + Duu + fu

Fig. 3

 3 Fig. 3 Neural oscillator scheme

3 Fig. 4

 34 Fig.[START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] Stable periodic steady-state oscillation of the neural oscillator projected into the (x 1 , x 3 ) plane

5 Fig. 5

 55 Fig. 5 Neural system. Error in the steady-state solution computation obtained with different discretization techniques and with N ∈ [50, 600]: maximum error in the x 2 trajectory

Fig. 8

 8 Fig. 8 Nyquist diagram of the transfer function of the dynamical system (continuous line) together with the representation of the negative reciprocal of the deadzone describing function (dashed line)

3 Fig. 9 5 Fig. 10

 39510 Fig.[START_REF] Li | A new method to determine the periodic orbit of a nonlinear dynamic systems and its period[END_REF] Unstable periodic oscillation for the deadzone feedback system: steady-state solution computed by using the MQCP approach (continuous line) and diverging time-stepping simulation (dotted line)

5 Fig. 11

 511 Fig. 11 Deadzone feedback system. Error in the period computation obtained with different discretization techniques and with N ∈ [100, 5400]

Fig. 13

 13 Fig.[START_REF] Leine | Stick-slip vibrations induced by alternate friction models[END_REF] Periodic steady-state solution of the stick-slip system. The graphs show the elongation of the spring, x 1 , the velocity of the mass, x 2 , and the friction force, λ. The solution is obtained by fixing N = 1000, θ = 0.5 and γ = 1. The corresponding value of the period is 4.9355 s

Fig. 14 Fig. 15

 1415 Fig. 14 Stick-slip system. Error in the steady-state solution computed with different discretization techniques and with N ∈ [90, 1080]: maximum error in the x 1 trajectory

Fig. 16 Fig. 17

 1617 Fig.[START_REF] Lu | A modified incremental harmonic balance method for rotary periodic motions[END_REF] Periodic steady-state solution of a repressilator network consisting of three genes. The graph shows the evolution of the variables x 1 , x 2 and x 3 which represent the concentration of proteins Lacl, Tetr and Cl, respectively. The solution is obtained by fixing N = 1080, θ = 0.5 and γ = 1. The corresponding value of the period is T = 15.1116 s

Fig. 18

 18 Fig.[START_REF] Brambilla | MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits[END_REF] Repressilator network. Error in the steady-state period computation obtained with different discretization techniques and with N ∈[START_REF] Sessa | Computing period and shape of oscillations in piecewise linear Lur'e systems: a complementarity approach[END_REF] 1080]

Fig. 19 -Fig. 20

 1920 Fig.[START_REF] Duarte | Fractional describing function of systems with Coulomb friction[END_REF] Relay feedback system: effect of a grazing-sliding bifurcation. Periodic steady-state trajectory of the variable x 1 : a for ζ = 0.082 it shows six pieces of sliding, c and for ζ = 0.091 the trajectory grazes the switching surface at time instants 0.3475

  cal system[START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] (continuous line) together with the representation of the negative reciprocal of the relay describing function (dashed line)

		0.2
	x 1	0
	-0.2
		0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8
						time [s]
		0.1
	x 2	0
	0 -0.1	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8
						time [s]
		0.02
	x 3	0
	-0.02
		0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8
						time [s]