Javier Rodríguez-Rodríguez

Almudena Casado-Chacón

Daniel Fuster

Physics of beer tapping

The popular bar prank known in colloquial English as beer tapping consists in hitting the top of a beer bottle with a solid object, usually another bottle, to trigger the foaming over of the former within a few seconds. Despite the trick being known for long time, to the best of our knowledge, the phenomenon still lacks scientific explanation. Although it seems natural to think that shock-induced cavitation enhances the diffusion of CO2 from the supersaturated bulk liquid into the bubbles by breaking them up, the subtle mechanism by which this happens remains unknown. Here we show that the overall foaming-over process can be divided into three stages where different physical phenomena take place in different time-scales, namely: bubble-collapse (or cavitation) stage, diffusion-driven stage and buoyancy-driven stage. In the bubble-collapse stage, the impact generates a train of expansion-compression waves in the liquid that leads to the fragmentation of pre-existing gas cavities. Upon bubble fragmentation, the sudden increase of the interface-area-to-volume ratio enhances mass transfer significantly, which makes the bubble volume grow by a large factor until CO2 is locally depleted. At that point buoyancy takes over, making the bubble clouds rise and eventually form buoyant vortex rings whose volume grows fast due to the feedback between the buoyancy-induced rising speed and the advection-enhanced CO2 transport from the bulk liquid to the bubble. The physics behind this explosive process sheds insight into the dynamics of geological phenomena such as limnic eruptions.

Understanding the formation of foam in a supersaturated carbonated liquid after an impact on the container involves a careful physical description of a number of processes of great interest in several areas of Physics and Chemistry. In order of appearance in this problem: propagation of strong pressure waves in bubbly liquids, bubble collapse and fragmentation, gas-liquid diffusive mass transfer and the dynamics of bubble-laden plumes and vortex rings. All these phenomena are observed, for instance, in the explosive formation of foam occurring in a beer bottle when it is tapped on its mouth, an effect known as beer tapping. In this Letter, we will use this effect as a convenient system to quantitatively describe the interaction between the processes mentioned above that ultimately leads to the explosive formation of foam that occurs in gas-driven eruptions [START_REF] Zhang | Dynamics of lake eruptions and possible ocean eruptions[END_REF]. As a consequence of the broad range of phenomena taking part in the overall process, the better understanding of the foam forming process in supersaturated liquids finds application in various fields of natural sciences and technology where similar gas-driven eruptions occur. The dynamics of limnic [START_REF] Zhang | Dynamics of lake eruptions and possible ocean eruptions[END_REF][START_REF] Zhang | Dynamics of CO2-driven lake eruptions[END_REF] or explosive volcanic [START_REF] Cashman | How volcanoes work: a 25 year perspective[END_REF][START_REF] Mangan | Gas evolution in eruptive conduits: combining insights from high temperature and pressure decompression experiments with steady-state flow modeling[END_REF] eruptions and the formation of flavour-releasing aerosols by bursting Champagne bubbles [5] are just a few examples. Mott and Woods [START_REF] Mott | A model of overturn of CO2 lakes triggered by bottom mixing[END_REF] have triggered a chain reaction in a stably-stratified tank containing a deep layer of CO 2 -saturated lemonade and a shallower layer of fresh water by spilling a gravity current of salt grains along the bottom of the tank. This example shows that the dynamics of CO 2 bubbles in daily-life liquids can be used to explain complex natural phenomena such as limnic eruptions.

To understand how the processes described above interact to lead to the foaming-up of beer, we have carried out an experimental investigation impacting commercial beer bottles under well-controlled repeatable conditions (see Supplementary Material). In particular, bubbles have been generated at a fixed location far from the bottle walls by focusing a laser pulse into the bulk liquid. In this way, we avoid the variability in the formation of bubbles upon the impact caused by the arbitrary distribution of nucleation sites [START_REF] Brennen | Cavitation and bubble dynamics[END_REF], thus ensuring that a bubble with a known initial size is always present in the measurement volume. By recording the evolution of these gas bubbles with a high-speed camera and the liquid pressure temporal evolution with a hydrophone we provide qualitative and quantitative analyses of the various processes that develop during the foam formation. Thus, we divide the overall foaming-over process into well-differentiated stages controlled by different physical mechanisms. More importantly, we show experimental evidence supporting the explanation given for each step of the outgassing process.

The chain of events that ultimately leads to the foaming-up of beer is triggered by a sudden impact on the top of the bottle, which generates a compression wave that propagates through the glass towards the bottom as predicted by the classical theory of impact on solids [START_REF] Goldsmith | Impact: The Theory and Physical Behaviour of Colliding Solids[END_REF]. When the wave reaches the base of the bottle, it is partially transmitted to the liquid as an expansion wave that travels towards the free surface, where it bounces back as a compression wave. A simple model for the impact problem [START_REF] Goldsmith | Impact: The Theory and Physical Behaviour of Colliding Solids[END_REF] implies that the stiffer the bottle, the more efficient is the transmission of the expansion wave to the liquid. Thus, the shock is more efficiently transmitted into the liquid in the case of a glass container than in a softer bottle (e.g. plastic), although the expansion wave is still generated, albeit its smaller amplitude. The train of waves transmitted to the liquid bounce back and forth several times until it damps out. Figure 1 shows snapshots of the first instants after the impact to illustrate the effects of the first expansion-compression cycle. It can be seen how bubbles start to expand first near the bottom whereas, at approximately t ≈ 124 µs, those located near the free surface begin to shrink. The train of rarefactioncompression waves drives the fragmentation of most of the existing gas pockets during the first wave cycles. Figure 2 illustrates a typical time evolution of the bubble radius measured in our experiments along with snapshots showing the bubble at different relevant instants of the expansion-collapse process (see also Movies S1 and S2 in Supplementary Material). This first stage of the overall foaming-up process lasts of the order of the acoustic time of the liquid volume, usually t ac = 2H/c ≈ 0.2 ms, assuming a typical liquid height H ≈ 10 cm and a speed of sound around c ≈ 1000 m/s. Notice that, since waves propagating in bubbly liquids are strongly damped, the intensity of successive rebounded waves decays rapidly thus they are less likely to cause bubble collapse.

Similarly to what happens in the generation of medical ultrasound contrast agents through sonication (ref. [START_REF] Feinstein | Two-dimensional contrast echocardiography. I. In vitro development and quantitative analysis of echo contrast agents[END_REF]) or, albeit in a more violent way, in sonoluminescence [START_REF] Brenner | Single-bubble sonoluminiscence[END_REF], it seems reasonable to attribute the break-up of the bubbles to a Rayleigh-Taylor instability [START_REF] Brenner | Single-bubble sonoluminiscence[END_REF]. The number of fragments, N , resulting upon the break up of a bubble cannot be measured due to the high void fraction of the resulting bubble cloud (Fig. 2d). Instead, an estimation of this number is obtained using the model of Brennen [START_REF] Brennen | Fission of collapsing cavitation bubbles[END_REF], based on the ideas put forward by several authors [START_REF] Prosperetti | Linear stability of a growing or collapsing bubble in a slightly viscous liquid[END_REF][START_REF] Ioss | Stability of a compressed gas bubble in a viscous fluid[END_REF]. Following this model (see Supplementary Material), the most unstable mode is given by

n m = 1 3 (7 + 3Γ m ) 1/2 -2 , (1)
with Γ m = ρR 2 R/σ, evaluated at the instant when the radius, R, is minimum, ρ the fluid density and σ the liquid-gas surface tension. The size of the fragments is expected to be of the order of R/n m , thus the number of fragments generated is N ≈ n 3 m . For the typical bubble sizes and pressure wave amplitudes used in these experiments, we find a most unstable mode of the order of n m ≈ 10 2 and a number of fragments N ≈ 10 6 .

As a consequence of the fast bubble collapse and breakup, right after the implosion the total gas-liquid interfacial area increases by a factor of the order of N 1/3 . This sudden increase of the interfacial area leads to a second stage where the clouds of bubble fragments grow rapidly as a result of the diffusion of carbonic gas into the newly created cavities. This stage can be modelled using the classical theory of bubble growth in supersaturated media (Ref. [START_REF] Epstein | On the stability of gas bubbles in liquid-gas solutions[END_REF]). Under the reasonable assumption that the cloud grows as the sum of its components, this theory states that the cloud size, L c , follows

L c = L 0 + αN 1/3 F ∆C ρ g κt π , (2)
where ∆C is the difference between the concentration of carbonic gas in the bulk liquid and the saturation value, κ its diffusivity, ρ g the density of the gas inside the bubbles, α a dimensionless constant and F (x) a known function (Supplementary Material). Taking the estimated number of fragments generated during the collapse of a single bubble, N ≈ 10 6 , we expect the radius of the bubble cloud to grow about 100 times faster than a single bubble with the same volume than the cloud. In fact, this magnitude represents an upper bound, since those fragments at the center of the cloud will grow more slowly, due to their limited access to CO 2 . Initially, the growth rate scales roughly as t 1/2 albeit exhibiting some oscillations caused by cycles of expansions and compressions that are not yet attenuated (see the stage labelled as "diffusiondriven" in Fig. 3, blue squares in the upper panel). This diffusion-driven stage ends when carbon dioxide is locally depleted and thus the cloud's size significantly moderates its growth. In the example of figure 3, this occurs at about t ≈ 10 ms (stage labelled as "depletion"). To avoid the noise introduced by the acoustic waves at short times we have performed additional experiments by focusing a highenergy laser pulse inside the liquid to trigger the formation of a bubble cloud by laser-induced cavitation [START_REF] Ohl | Luminescence from spherically and aspherically collapsing laser induced bubbles[END_REF] (see Supplementary Material). In the focal region, the laser generates a dense bubble cloud that initially grows as the square root of time, which is consistent with a purely diffusive growth (Fig. 3, red squares in the upper panel).

The rapidly growing bubble clusters act as buoyancy sources that lead to the formation of bubble-laden buoyant vortex rings in time scales of order t g ∼ (L/g) 1/2 (Fig. 3c), very much like a localized release of heat forms a thermal [START_REF] Turner | Buoyant vortex rings[END_REF]. As the vortices rise through the liquid, the advection due to their self-induced velocity and the mixing caused by their vortical motion contribute to enhance the transport of CO 2 to the bubbles. In turn, this results in a growth-rate faster than that found for pure diffusion, namely t 1/2 . Indeed, the cloud's size grows roughly as t 2 during this stage (Fig. 3, upper plot). Moreover, as a consequence of the continuous generation of gas volume inside the vortex, the velocity approaches a constant value (Fig. 3, lower plot) instead of decreasing as happens in buoyant vortex rings induced by the release of a fixed amount of buoyancy. For instance, in thermals originated by the sudden release of a fixed amount of heat, the velocity decays as t -1/2 due to the entrainment of colder fluid [START_REF] Turner | Buoyant vortex rings[END_REF][START_REF] Morton | Weak thermal vortex rings[END_REF].

Conversely, in the bubble-laden plumes studied here, the feedback between buoyancy-driven rising motion and gas-volume generation results in a nearly constant speed. This behavior is similar to that found in the the socalled autocatalytic vortex rings or plumes [START_REF] Rogers | Buoyant Plumes and Vortex Rings in an Autocatalytic Chemical Reaction[END_REF] where buoyancy is continuously produced by a chemical reaction that yields products less dense than the reactants. These buoyancy-driven chemically-reacting flows appear, for instance, in the combustion of flame balls in microgravity conditions. Interestingly, they are also relevant in some explosion scenarios for Type Ia supernovae [START_REF] Vladimirova | Model flames in the Boussinesq limit: Rising bubbles[END_REF][START_REF] Rogers | The heads and tails of buoyant autocatalytic balls[END_REF].

The analogy between the bubble-laden plumes observed here and the autocatalytic plumes described in the literature extends also to their morphology. Panel c in Fig. 3 shows one of these bubble-laden plumes when it is well developed (see also Movie S3). The plume consists of a vortex with a nearly spherical cap with a thin conduit that ascends more slowly, features observed in the plumes driven by autocatalytic chemical reactions [START_REF] Rogers | Buoyant Plumes and Vortex Rings in an Autocatalytic Chemical Reaction[END_REF][START_REF] Rogers | The heads and tails of buoyant autocatalytic balls[END_REF].

It should be pointed out that, among all the stages of the foaming-up process, this one is the most effective in terms of the amount of liquid outgassing as a result of its self-accelerating nature. This stage starts at times of the order of tens of milliseconds and concludes when the plumes reach the size of the liquid volume, usually of the order of a second.

Remarkably, the behavior of the bubble-laden vortex rings during the diffusion-driven and buoyancy-driven stages is independent of the mechanism used to generate the initial bubble cloud. Figure 3 depicts the evolution of the bubble cluster size and velocity of a bubble cluster originated by laser-induced cavitation. The size and velocity follow the same scaling laws as the vortex created by the pressure-induced bubble implosion. This suggests that similar explosive CO 2 outgassing processes driven by the formation of these bubbly plumes may be initiated by other physical mechanisms generating dense bubble clouds such as the introduction of new bubble nucleation sites [START_REF] Coffey | Diet coke and mentos: what is really behind this physical reaction?[END_REF][START_REF] Woods | Turbulent Plumes in Nature[END_REF] or a sudden change on the saturation conditions occurring either globally [START_REF] Mangan | Gas evolution in eruptive conduits: combining insights from high temperature and pressure decompression experiments with steady-state flow modeling[END_REF] or locally [START_REF] Zhang | Dynamics of CO2-driven lake eruptions[END_REF][START_REF] Woods | Turbulent Plumes in Nature[END_REF]. In fact, our observations suggest that, once the plume is initiated, its dynamics do not seem to depend on the particular initiation mechanism. Thus, one of the main conclusions of this study is that the dynamics of these bubble-laden self-accelerating plumes moving in supersaturated media may partly explain the explosive behavior of systems such as limnic and explosive volcanic eruptions where current models typically neglect the role of these autocatalytic structures [START_REF] Zhang | Dynamics of lake eruptions and possible ocean eruptions[END_REF][START_REF] Cashman | How volcanoes work: a 25 year perspective[END_REF][START_REF] Mangan | Gas evolution in eruptive conduits: combining insights from high temperature and pressure decompression experiments with steady-state flow modeling[END_REF][START_REF] Mott | A model of overturn of CO2 lakes triggered by bottom mixing[END_REF][START_REF] Woods | Turbulent Plumes in Nature[END_REF].

Finally, two side effects induced by the development of the bubbly plumes must be mentioned here attending to their relevance in the global degassing process in the case of the bottle. Firstly, due to the finite size of the container, a global recirculating motion is generated that drags bubbles from near the free surface deep into the bulk liquid, thus increasing their residence time in the flow and allowing them to grow for longer times. Secondly, the flow induced inside the bottle speeds-up also the growth of gas cavities attached to the walls [START_REF] Jones | Bubble nucleation from gas cavities -a review[END_REF][START_REF] Liger-Belair | Modeling nonclassical heterogeneous bubble nucleation from cellulose fibers: application to bubbling in carbonated beverages[END_REF] through the enhancement in the transport of carbon dioxide towards these cavities, that otherwise would only grow by diffusion. A very similar effect is probably behind the long time scales involved in limnic eruptions. In these phenomena, bubbly plumes form that keep entraining CO 2 -saturated water from the bottom of the lake until it is almost depleted of this gas due to the global overturning flow that they induce in the lake [START_REF] Mott | A model of overturn of CO2 lakes triggered by bottom mixing[END_REF]. Altogether, the chain of effects described in this Letter leads to the fast appearance of foam that has granted beer tapping its popularity.

We would like to thank Dr. Maylis Landeau for her useful comments on the potential applications of the present study. The authors are indebted to Profs. Norman Riley, José Manuel Gordillo, Alejandro Sevilla and Carlos Martínez-Bazán for their insightful comments. We acknowledge the support of the Spanish Ministry of Economy and Competitiveness through grant DPI2011-28356-C03-02.

SUPPLEMENTARY MATERIAL

Experimental techniques:

In a first set of experiments, the pressure fluctuations induced in the liquid upon the impact were recorded with a hydrophone. In these experiments, the bottle was filled with deionized water, to delay the formation of cavitation bubbles at the hydrophone's surface. In a second set of experiments, actual beer was used (CO 2 concentration: 5.0±0.4 g/l, measurement provided by the manufacturer). The evolution of bubbles existing in the liquid previously to the impact was recorded with a high-speed camera. To ensure that a bubble is present at the camera's focal point, bubbles were generated by focusing a low-energy YaG laser pulse inside the bulk liquid with a convergent lens before the impact. In this way, we avoid the variability in the location and size of the bubble formed upon the impact due to the randomness of the distribution of nucleation sites found in all liquids [START_REF] Brennen | Cavitation and bubble dynamics[END_REF]. In these two experimental sets, the bottle was impacted by dropping a brass weight (90 g.) in a controlled and repeatable manner from a height of 25 mm. (see Fig. 4). Moreover, the bottle is held by the neck using a rigid clamp that prevents the net motion of the bottle after the impact and keeps the bottom above the ground, thus preventing the partial transmission of the wave to the bench. The bottle is held in place with a clamp around its neck that prevents its motion after the impact. To ensure a good repeatability, the weight falls along a rail that guarantees that its lower side remains horizontal during the falling. The YaG laser pulse used to generate both the individual bubbles (low-energy pulse) and the bubble clouds (high-energy pulse) is also depicted.

Rail

Finally, a third kind of experiment was performed where, instead of hitting the bottle mechanically, a high-energy laser pulse was focused inside the beer to induce cavitation, similarly as described in [START_REF] Ohl | Luminescence from spherically and aspherically collapsing laser induced bubbles[END_REF], thus creating a bubble cloud that serves as a seed to eventually trigger the formation of a bubble-laden plume. In all the sets, the projected area of the bubble cluster, A, was measured using custom-made software. The cluster size, L c , was then calculated as the diameter of a circle with the same area, namely L c = 4A/π.

Estimation of the number of bubble fragments: Upon bubble collapse, the high void fraction of the resulting cloud of fragments precludes the application of any measurement technique to directly determine the number of those fragments. To get at least an order of magnitude estimation of that number, we apply the model proposed by Brennen [START_REF] Brennen | Fission of collapsing cavitation bubbles[END_REF]. Briefly, this model uses the equation for the amplitude, a(t), of a spherical harmonic perturbation of order n > 1,

ä+ 3 R Ṙ ȧ-(n -1) R R -(n -1)(n + 1)(n + 2) σ ρR 3 a = 0, (3)
where ρ is the liquid density and σ the liquid-gas surface tension. Notice that, to determine the evolution of the perturbations, the time-history of the bubble radius, R(t), is needed. As a further simplification, since the growth of the surface instabilities is very fast compared to the bubble radial dynamics, both the radius, R, and its acceleration, R, are considered constant during the final instants of the break-up process. Thus, it is now possible to calculate the fastest-growing mode, n m , i.e. that maximizes the expression in brackets, [...], in equation 3:

n m = 1 3 (7 + 3Γ m ) 1/2 -2 , (4)
where Γ m = ρR 2 R/σ, evaluated using the minimum bubble radius, R min and the acceleration at the time at which that radius is reached. Consequently, Brennen's model allows us to compute the fastest-growing mode using only values of the radius and its acceleration at the collapse time.

Nonetheless, the short duration of the collapse stage impedes recording the process with enough time resolution to estimate either the minimum bubble radius or the acceleration during this stage (see Figs. 2c and2d). Therefore, to obtain these parameters to apply Brennen's model [START_REF] Brennen | Fission of collapsing cavitation bubbles[END_REF], the evolution of the bubble radius has been computed numerically by integrating the Rayleigh-Plesset equation:

ρR R + 3 2 ρ Ṙ2 -p 0 + 2σ R 0 R 0 R 3γ + +p 0 + 2σ R + 4µ Ṙ R = -p w (t), (5)
where µ = 10 -3 kg/(m•s) is the liquid viscosity, γ = 1.304 is the heat capacity ratio of the gas (CO 2 ) and p w (t) is the pressure pulse induced in the system, given by p w (t) = -p A sin (2πt/T ), with amplitude p A = 100 kPa and period T = 0.24 ms. These values were measured with a hydrophone at the bubble's location. Moreover, the following values have been used: the CO 2 -water surface tension, σ = 0.0434 N/m, the water density, ρ = 10 3 kg/m 3 , and the ambient pressure, p 0 = 10 5 Pa.

Figure 2 shows the good agreement between the numerical and the experimental results before the break up, justifying the usage of the numerical parameters to estimate the number of fragments resulting from it. For instance, from the experiment shown in figure 2, we obtain n m ≈ 144. This value ranges between 120-160 in all our experiments, with a maximum different between the maximum radius calculated numerically and the measurement always within 10% or less.

Modelling the growth rate of a bubble cloud: To get an estimation of the growth rate of the bubble cloud due to the diffusion of CO 2 before buoyancy becomes important (diffusion-stage), we assume that the cluster's volume grows as the sum of the volume of its components. The time evolution of the radius of a fragment, R f (t), is given be the Epstein-Plesset equation (Ref. [START_REF] Epstein | On the stability of gas bubbles in liquid-gas solutions[END_REF]),

dR f dt = κ∆C ρ g 1 R f + 1 √ πκt , (6)
where κ is the diffusivity of CO 2 in water, ρ g the density of the gas inside the bubble and ∆C the difference between the concentration of CO 2 far away from the cloud and that at the bubble's surface. Neglecting the capillary pressure, the solution to the equation is well approximated by

R f (t) ≈ R f,0 + F ∆C ρ g κt π , (7)
with R f,0 the initial radius and F (x) =

x 1 + 1 + 2π/x (ref. [START_REF] Epstein | On the stability of gas bubbles in liquid-gas solutions[END_REF]). Notice that the second term in this expression corresponds to the solution to Eq. ( 7) at long times. Here we have added the first term, the initial radius, to obtain an expression that yields small relative errors for both long and short times. Assuming that the N fragments are equally sized, the volume of the cloud can be written as

L 3 c = α 3 N R f,0 + F ∆C ρ g κt π 3 , (8)
where the dimensionless constant α accounts for the fact that the cluster contains some amount of liquid separating the bubbles. Taking the cube root of the expression, and denoting the initial cluster size by L 0 = αN 1/3 R f,0 , we get

L c = L 0 + αN 1/3 F ∆C ρ g κt π . (9)
Finally, we would like to comment on the applicability of the Epstein-Plesset equation to this problem. This equation describes the growth of an isolated bubble, whereas in this case bubbles belong to a dense cluster. Therefore we expect that, at some point, the CO 2 must be depleted, with only those bubbles closer to the clouds' surface being still able to grow. This is consistent with the behavior observed in figure 3b, where the growth rate of the cluster decreases notably at about 10 ms.

MOVIE LABELS:

Movie S1. Bubble fragmentation 1: fragmentation of a bubble as a result of the passage of the train of expansion-compression waves. This is the same bubble shown in Fig. 2. Scale bar: 500 µm. Movie S2. Bubble fragmentation 2: the fragmentation of a gas bubble is shown, together with a plot of its instantaneous radius. Notice how, after the second compression, the bubble breaks-up and forms a cloud of bubble fragments.

Movie S3.

Bubble plume development: this movie shows the evolution of a bubble cloud during the buoyancy-driven stage. Notice the large increase of the foam volume. This plume corresponds to the blue symbols of Fig. 3. Scale bar: 1 mm. Movie S4. Global view of the bottle after the impact: the movie shows the development and growth of the bubble plumes after hitting the bottle. To make more evident the process, the formation of bubbles has been forced by placing a metal disk at the bottom of the bottle, what effectively introduces a large number of nucleation sites at the bottom.

FIG. 1 .

 1 FIG. 1. Sequence of images corresponding to the instants right after the pressure wave starts to propagate through the beer. To ensure a continuous bubble cloud, a metal disk was placed at the bottom of the bottle in this experiment, what effectively introduces a large number of nucleation sites. Scale bar: 10 mm.

FIG. 2

 2 FIG.2. A bubble of initial radius R0 ≈ 180 µm, induced by a laser pulse nearly a second before the impact, grows after the passage of the expansion wave reaching the bubble location at instant (a). After the instant of maximum radius (b), the bubble collapses at some point between frame (c) and the previous one, turning into a bubble cloud. Panel (d) shows the bubble cloud 0.1 ms after instant (c). The Rayleigh-Plesset equation has been integrated numerically (red solid line, see Supplementary Material) for a bubble subjected to a pressure pulse with an amplitude, pA = 100 kPa, and period, T = 0.24 ms, measured at the bubble's location with a hydrophone. Notice that, after the implosion, t ≈ 0.22 ms, the Rayleigh-Plesset equation no longer describes the behavior of the bubble since it is only valid for a single bubble, not for a bubble cloud. (see Movies S1 and S2).

FIG. 3 .

 3 FIG.3. Upper plot: Time evolution of the size of a bubble cloud, Lc, after the shock-induced collapse of an already existing bubble (blue) and after a laser-induced bubble implosion (red). Three different stages have been marked qualitatively: diffusion driven, depletion and buoyancy-driven. Black solid lines depict the scaling laws Lc ∼ t 1/2 (diffusion-driven) and Lc ∼ t 2 (buoyancy-driven). Lower plot: Time evolution of the top of the bubble cloud corresponding to the same cases of the upper plot. Notice how, at long times, the plume approaches a steady rising velocity. Letters denote the instants corresponding to the images in the lower pannels (scale bar: 1 mm). (See movie S3)

FIG. 4 .

 4 FIG.4. Sketch of the experimental setup. The bottle is held in place with a clamp around its neck that prevents its motion after the impact. To ensure a good repeatability, the weight falls along a rail that guarantees that its lower side remains horizontal during the falling. The YaG laser pulse used to generate both the individual bubbles (low-energy pulse) and the bubble clouds (high-energy pulse) is also depicted.


Mara Virginia Galmarini

Anne-Laure Loiseau

Michel Visalli

Pascal Schlich

How about some cheese

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

Cheese intake

 Pacherenc: high score ratings did not change, in spite of the impact of cheese on the duration of dominance of 4 descriptors.

 Sancerre: all cheeses reduced the duration of perceived astringency, but only Comté increased the liking scores.

 Bourgogne: all cheeses reduced the dominance of astringency and some had an impact on red fruits aroma. How about some cheese with that wine? Use of Multi-Intake Temporal Dominance of Sensations to evaluate the influence of cheese on wine perception and appreciation. The effect of cheese on wine perception was evaluated by comparing the TDS profile of the 3 rd sip of wine with no cheese, to those obtained after cheese intake.

Details on data analysis

Multi-intake Temporal Dominance of Sensations (TDS)

alternated with an hedonic test.

Wines were evaluated by multi-intake TDS according to the protocol described in the figure below. For each sip of wine the TDS profile was done; and then an hedonic score was given. Before taking a second sip of wine, consumers eat a piece of cheese. Right after finishing eating they evaluated a second sip. The cycle was repeated 3 times. The 16 wine-cheese combinations were evaluated over 4 sessions. As a reference, the same task (TDS + hedonic test) was done for each wine without cheese intake.

Some "wine-cheese couples" have a good reputation while others are discouraged of being consumed together, based on terroir or tradition. However, there is no proof on how one cheese could influence the sensory perception of a wine.

Hedonic evaluation CONCLUSION

 Cheese had an impact (product dependent) on the temporal perception and appreciation of wine, observed thanks to multi-intake TDS.  As a general rule, cheese reduced the duration of dominance of astringency in wine.

 Wine appreciation was improved by the presence of cheese; in other words, some wines should rather be drunk with cheese than alone. Preference was not reduced by any of the cheeses.

Wines and Cheeses: chosen on usual gastronomic recommendations

Protocol for wine evaluation together with cheese intake

MATERIALS & METHODS

RESULTS -Focus on Madiran wine

Reduction of duration of dominance of sour and astringency after eating cheese was significant and equal for all cheeses. The opposite was observed for red fruits which increased with most cheeses.

Bandplots by attribute for wine with no cheese and after Comté, Crottin de Chavignol, Epoisses and Roquefort.

Differences on sequentiality of significant dominant sensations.

TDS evaluation of one sip 3 sips

The same descriptors appear as dominant in most cases except for: alcoholic (dominant only after cheese consumption)

and sour (dominant when no cheese was consumed and after Roquefort). There were also some changes on sequentiality E.g.: woody was dominant in the middle and the end of the sip for Madiran but at the beginning for Madiran+roquefort.

Differences on duration of dominance 1,2 Differences on appreciation 3

When drunk alone, wine appreciation decreased from Sip 1 to 3. An opposite behavior was observed after cheese intake, regardless of the type of cheese.

Summary for the other wine-cheese combinations

Significant changes in duration of dominance

The panel: consisted of 31 frequent wine and cheese consumers somewhat familiar with TDS.

Crottin de Chavignol Goat cheese

Roquefort Blue cheese

Comté

1

  -Representation of Canonical Variate analysis (CVA) of durations. 2-ANOVA/MANOVA ; Duration = Product +Subject 3-Done by wine, comparing the effect of cheese intake, using sip as repeated measure. ProcMixed, SAS ®

  5; **p<0.01; ***p<0.001. Different letters shows significant differences according to an LSD test.(.)p<0.1; *p<0.5; **p<0.01; ***p<0.001. Different letters shows significant differences according to an LSD test.

1,2 Galmarini, Mara V.*; 1 Loiseau, Anne-Laure, 1 Visalli, Michel; 1 Schlich, Pascal 1

  Centre des Sciences du Goût et de l'Alimentation (CSGA), INRA, Dijon, France 2 Consejo Nacional de Investigaciones Científico Tecnológicas (CONICET), Argentina


Antoine Liutkus
email: antoine.liutkus@inria.fr

Kazuyoshi Yoshii
email: yoshii@kuis.kyoto-u.ac.jp

A DIAGONAL PLUS LOW-RANK COVARIANCE MODEL FOR COMPUTATIONALLY EFFICIENT SOURCE SEPARATION

Keywords: Blind source separation, nonnegative matrix factorization, positive semidefinite tensor factorization, low-rank approximation

This paper presents an accelerated version of positive semidefinite tensor factorization (PSDTF) for blind source separation. PSDTF works better than nonnegative matrix factorization (NMF) by dropping the arguable assumption that audio signals can be whitened in the frequency domain by using short-term Fourier transform (STFT). Indeed, this assumption only holds true in an ideal situation where each frame is infinitely long and the target signal is completely stationary in each frame. PSDTF thus deals with full covariance matrices over frequency bins instead of forcing them to be diagonal as in NMF. Although PSDTF significantly outperforms NMF in terms of separation performance, it suffers from a heavy computational cost due to the repeated inversion of big covariance matrices. To solve this problem, we propose an intermediate model based on diagonal plus low-rank covariance matrices and derive the expectation-maximization (EM) algorithm for efficiently updating the parameters of PSDTF. Experimental results showed that our method can dramatically reduce the complexity of PSDTF by several orders of magnitude without a significant decrease in separation performance.

INTRODUCTION

A major approach to blind source separation of single-channel audio signals in the last decade is to use nonnegative matrix factorization (NMF) and Wiener filtering by assuming that all the time-frequency bins in the short-time Fourier transform (STFT) domain are independent from each other [START_REF] Benaroya | Audio source separation with a single sensor[END_REF][START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF][START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. In NMF, given a set of nonnegative vectors as input data, each vector is approximated by the weighted sum of nonnegative basis vectors. In audio source separation, a mixture spectrum (magnitude, power, or α-fractional power spectrum [START_REF] Liutkus | Generalized Wiener filtering with fractional power spectrograms[END_REF][START_REF] Liutkus | Alpha-stable matrix factorization[END_REF]) is Each nonnegative vector (power spectrum) is approximated by the weighted sum of fewer nonnegative vectors Each positive semidefinite matrix (complex covariance matrix) is approximated by the weighted sum of fewer positive semidefinite matrices : complex spectrum Fig. 1. Comparion between LD-PSDTF and IS-NMF.

approximated as the weighted sum of source spectra at each frame. Many variants of NMF have been proposed by designing various cost functions that evaluate the approximation error at each time-frequency bin [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-way Data Analysis and Blind Source Separation[END_REF]. Among these, NMF based on the Kullback-Leibler (KL) [START_REF] Smaragdis | Static and dynamic source separation using nonnegative factorizations: A unified view[END_REF] or β-divergence [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF][START_REF] Févotte | Algorithms for nonnegative matrix factorization with the beta-divergence[END_REF] have been empirically known to work well. NMF based on the Itakura-Saito (IS) divergence is known to be theoretically justified for audio source separation under the assumption that all the coefficients of the sources STFT are independent and isotropic Gaussian [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF]. In this context, it can indeed be interpreted as the maximum likelihood estimation of the sources variances, that are called power spectral densities. Separation through posterior inference of the Gaussian sources conditionally on the mixture then coincides with a frame-wise Wiener filter. A serious limitation of this approach is that, while the magnitude of the input mixture STFT is partitioned to the different sources, all sources get the same original phase information. It is thus difficult to resynthesize high-quality time-domain source signals.

To overcome this limitation, positive semidefinite tensor factorization (PSDTF) based on the log-det (LD) divergence has recently been proposed [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF] by assuming that the source complex spectrum at each frame is multivariate complex Gaussian distributed. In PSDTF, given a set of positive semidefinite (PSD) matrices as input data, each PSD matrix is approximated by the weighted sum of PSD basis matrices (Fig. 1). Focusing on the fact that the positive semidefiniteness of matrices is an extended concept of the nonnegativity of vectors, LD-PSDTF is a natural extension of IS-NMF. More specifically, each PSD matrix of input data is given by calculating the product of the complex spectrum and its conjugate transpose. NMF, on the other hand, focuses on only the diagonal elements of the PSD matrix, i.e., nonnegative vector (power spectrum), by ignoring the correlations between frequency bins. PSDTF can estimate the complex spectrum of each source with phase information by processing jointly all the frequencies in an interdependent manner. Since such phaseaware frequency-domain decomposition was shown to correspond to time-domain decomposition, PSDTF can recover high-quality time-domain source signals. Note that recentlyproposed phase-aware or time-domain decomposition methods based on additivity of complex spectra or audio signals do not take into account inter-frequency correlations [START_REF] Badeau | Gaussian modeling of mixtures of nonstationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Févotte | Low-rank time-frequency synthesis[END_REF][START_REF] Kameoka | Multi-resolution signal decomposition with time-domain spectrogram factorization[END_REF] (high-resolution NMF [START_REF] Badeau | Gaussian modeling of mixtures of nonstationary signals in the time-frequency domain (HR-NMF)[END_REF] deals with inter-frame dependency based on autoregressive modeling).

A critical problem of LD-PSDTF, however, lies in extremely large computational cost, which hinders its use in practice. When a mixture spectrogram with T frames and F frequency bins is analyzed, the time complexity of LD-PSDTF with K basis matrices is O(KT F 3 ) while that of IS-NMF with K basis vectors is only O(KT F ). LD-PSDTF is performed with a convergence-guaranteed iterative optimization method that needs to repeatedly calculate the inversion of very large matrices of size F × F . The matrix inversion costs O(F 3 ) and is numerically unstable because F (window size) is usually several thousands. If all the input and basis matrices are restricted to diagonal matrices by discarding the correlations between frequency bins, i.e., LD-PSDTF reduces to IS-NMF, the matrix inversion costs only O(F ) but the separation performance is deteriorated accordingly.

To solve this problem, we propose a constrained version of LD-PSDTF that represents each basis matrix as the sum of a diagonal matrix and a low-rank matrix (Fig. 2). This contribution is inspired by the factor analysis model presented in [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF], but goes further by adopting such a structure for each source. If the rank of the low-rank matrix is N F , the time complexity of our model is O(KT F 2 N ). If a basis matrix represents a harmonic musical instrument sound, the frequency bins corresponding to the harmonic partials are highly correlated to each other. This implies that the rank N can be reduced to around the number of harmonic partials without sacrificing separation performance too much. As either an expectation-maximization (EM) algorithm or an auxiliaryfunction-based (minorization-maximization) method can be used for IS-NMF, in this paper we derive an EM-based method with fast inversion of diagonal plus low-rank basis matrices. The main contribution of this paper is to propose an efficient and accurate approximation to full LD-PSDTF by leveraging the characteristics of basis matrices. This appears important because PSDTF has a great potential for fundamentally raising the performance of any audio analysis methods using NMF if the problem of computational cost is solved.

L k is F × N with N F . [v] is a diagonal matrix with diagonal v.

LOG-DET POSITIVE SEMIDEFINITE TENSOR FACTORIZATION (LD-PSDTF)

In this section, we review log-det positive semidefinite tensor factorization (LD-PSDTF) for audio source separation from the viewpoint of probabilistic modeling and newly derive the expectation-maximization (EM) algorithm for PSDTF.

Model formulation

All the signals are represented in the STFT domain, with F non-redundant frequency bins and T frames. In that domain, for one given frame t, the K source signals are F ×1 complex vectors written x kt . We take them as all independent and distributed according to the LD-PSDTF model as follows:

x kt ∼ N c (0, Y kt ) , (1)
where N c indicates the multivariate complex Gaussian distribution and Y kt is the covariance matrix. Y kt is assumed to the time-varying scaled version of a time-invariant F × F covariance template V k 0 as follows:

Y kt = h kt V k ,
with h kt ≥ 0 being an activation gain for source k at frame t.

As demonstrated in [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF], LD-PSDTF generalizes IS-NMF [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF] for which all entries of x kt are assumed to be independent, resulting in diagonal V k . The major advantage of LD-PSDTF is to overcome the limitation of this independence assumption. Since the perfect stationarity of the waveforms does not hold true and finite-size frames are used in practice, interfrequency covariances cannot be avoided. Now, we model the F × 1 mixture x t at frame t as the sum of the sources as follows:

x t = K k=1 x kt .
As a sum of independent Gaussian random vectors, x t also has a Gaussian distribution as follows:

x t ∼ N c 0, Y t = K k=1 Y kt .
This is the likelihood function of LD-PSDTF. Its maximization is known to be equivalent to the minimization of the log-det divergence between X t = x t x H t and Y t [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF]. Since a complex spectrum x t is the linear transform of a time-domain signal, such frequency-domain Gaussian modeling can be represented in the time domain.

Parameter estimation

To estimate the covariance templates V k as well as the activation gains h kt , we propose an generalized EM algorithm whose E-step and M-step are iterated until convergence. Note that we could use the auxiliary function approach presented in [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF]. Let Θ denote the set of all parameters. The whole procedure is summarized in Algorithm 1 and detailed below.

E-step

The posterior distribution of the sources x kt is computed. The posterior of Gaussian random variables is given by

x kt | Θ, x ∼ N c ( xkt , C kt ) , (2)
where the posterior mean xkt and covariance C kt are of dimension F × 1 and F × F , respectively, and given by

xkt = W kt x t , (3)
C kt = Y kt -W kt Y kt , (4)
with W kt being the F ×F Wiener gain for source k at frame t.

In full generality, it is given by

W kt = Y kt Y -1 t . (5)
The total posterior covariance of x kt is thus

E [x kt x kt | Θ, x] Σ kt = xkt x kt + C kt . (6)

M-step

Given the total posterior covariances Σ kt for the sources computed in the E-step, h kt and V k are alternatively updated. A straightforward option to estimate the activations for source k is to exploit the posterior distribution of x kt | Θ, x as achieved in (4). We then estimate h kt as:

h kt ← tr V -1 k Σ kt F . (7)
Note that this requires the inversion of the F × F matrix V k . This inversion needs to be done once per iteration though, and not for all frames. The update for the covariance template V k can also be derived when the posterior distribution for the sources has been computed as follows:

V k ← 1 T T t=1 Σ kt h kt . ( 8
)
Note that the likelihood function is not directly maximized in the M-step. Using the interdependent steps ( 7) and ( 8), the likelihood function is incrementally maximized.

Computational complexity

We now roughly estimate the time complexity of the EM algorithm. The iterative algorithm proposed here comprises sev-Algorithm 1 EM algorithm for LD-PSDTF. Input: x, K, initializations for V k and h kt Iterate until convergence:

1. E-step: compute all Σ kt as in ( 6) 2. M-step: update h kt with (7) and V k with (8). Separate: compute the sources xkt as in [START_REF] Liutkus | Generalized Wiener filtering with fractional power spectrograms[END_REF] eral demanding operations, in particular:

• Inversion of all F ×F covariance matrices Y t : O T F 3 .

• Computation of all W kt and W kt Y kt : O KT F 3 . Accounting only for these most demanding operations, the total computational complexity C b for each iteration of this baseline algorithm is thus 1 :

C b = O (K + 1) T F 3 . (9)
In this paper, we propose a refinement of LD-PSDTF that permits to dramatically reduce this computational load.

FAST LD-PSDTF

This section explain the proposed refinement of LD-PSDTF based on diagonal plus low-rank covariance matrices.

Diagonal plus low-rank approximation

Instead of leaving the covariance templates V k totally unconstrained as done in [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF], we assume that they can be approximated as the sum of low-rank and diagonal matrices (Fig. 2):

V k = [P k ] + L k [S k ] L k , (10)
where we write [v] as the diagonal matrix with vector v as its diagonal. L k is a F × N matrix, while P k and S k are F × 1 and N × 1 vectors, respectively, with N F . Note that we introduced S k to make parameter estimation easy (see Section 3.2.2).

We call this model (10) structured. It is for instance already reviewed in [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] under the name of factor analysis. Its rationale is to allow for the dependency structure between frequencies to be correctly explained by a limited amount of correlations. In other words, while P k may be roughly understood as the power spectral density (PSD) of the stochastic part of source k as in IS-NMF, L k [S k ] L k rather stands for its sinusoidal part, which introduces deterministic relations between the entries of x kt . Then, these two components are coupled because they share the activation gains h kt . Interestingly, if we take N = 0, this model reduces to IS-NMF and coincides with PSDTF for N = F . Now, taking [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF] as the covariance templates, the mixture covariance matrices Y t write:

Y t = K k=1 h kt P k + K k=1 h kt L k [S k ] L k .
The purpose of this paper is to study the consequences of model ( 10) over both computational load and performance in audio separation. For notational convenience below, and for a given p = 1, . . . , K, we define the F × F matrix M p,t as:

M p,t = K k=1 h kt P k + p k=1 h kt L k [S k ] L k .
With the structured model, the parameters P k , S k , and L k replace the unconstrained V k , leading to K (F + N + F N ) parameters for the templates instead of KF (F + 1) /2. The special structure [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF] chosen for the covariance templates V k may be understood as stating each source itself is the sum of a stochastic component x

kt = L k [h kt S k ] L k .
For the reestimation of these parameters, we modify the EM algorithm in the following way.

Paramter estimation

We explain an acceralated version of the generalized EM algorithm proposed in Section 2.2.

Modified E-step

We need to compute the total posterior covariances of x

Σ kt = Σ (s) kt + Σ (d) kt + Σ (ds) kt + Σ (ds) kt , (11)
where

Σ (ds) kt = x(d) kt x(s) kt -W (d) kt Y (s)
kt .

Modified M-step

The updates for P k are given by:

P k ← 1 T T t=1 diagΣ (s) kt h kt . ( 12
)
The updates for L k and S k are obtained straightforwardly by the truncated eigenvalue decomposition of the weighted average of Σ (d) kt as follows:

L k , S k ← eig 1 T T t=1 Σ (d) kt h kt , N , (13)
where eig (M, N ) gives the N -truncated eigenvalue decomposition of matrix M , i.e., the first N eigenvectors along with Algorithm 2 Inversion of the mixture covariance matrix Y t under the structured template model [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF]. 16) and [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF].

• M -1 0,t ← 1 k h kt P k • For k = 1 . . . K: Calculate (
• Return Y -1 t = M -1 K,t
their eigenvalues. Efficient randomized algorithms [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] may be used to compute this decomposition with a time complexity scaling as O F 2 N or even less.

Concerning the update for h kt , it is identical to [START_REF] Smaragdis | Static and dynamic source separation using nonnegative factorizations: A unified view[END_REF], except that we use the more efficient expression [START_REF] Badeau | Gaussian modeling of mixtures of nonstationary signals in the time-frequency domain (HR-NMF)[END_REF] for Σ kt .

Efficient matrix inversion

For inference in the EM algorithm (and also in the auxiliary function approach of [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF]), we need to compute the inverses of the F × F mix covariance matrices Y t . This step is one of the computational bottlenecks of the method. In this paper, we use the proposed structured model [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF] in conjunction with the Woodbury matrix identity [START_REF] Woodbury | Inverting modified matrices[END_REF], which states

(A + U CV ) -1 = A -1 -A -1 U C -1 + V A -1 U -1 V A -1 , (14)
for any matrices A, U , C, and V of appropriate size. The strategy we propose is to proceed iteratively.

Assuming we now have the inverse of M k-1,t , (which is a diagonal matrix for k = 1), we compute

M -1 k,t = (M k-1,t + h kt L k [S k ] L k ) -1 . (15)
Using Woodbury matrix identity ( 14) on [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF], we have

M -1 k,t = M -1 k-1,t -M -1 k-1,t L k Ω -1 k,t L k M -1 k-1,t , (16)
Ω k,t = 1 h kt S k + L k M -1 k-1,t L k , (17)
where 1 a denotes entry-wise inversion for a vector a. The procedure is iterated up to k = K, as shown in the Algorithm 2.

Note that the same identity ( 14) may be applied to compute V -1 k in a computationally efficient manner for the update of h kt in (7) as follows:

V -1 k = 1 P k (18)
I -L k 1 S k + L k 1 P k L k -1 L k 1 W k . (19)

Computational complexity

It can be shown that the computational complexity of each of the K iterations of the procedure presented in the algorithm box 2 in this structured case is

O N F 2 + F N + N 2 .
Then, considering our modifications to the generalized EM algorithm, we also need the add the following computationally most demanding operations: Keeping only the computations that need to be done for each frame, this brings the complexity of each iteration roughly down to

C s = O 3KT F 2 (N + 1) . (20)
In Fig. 3, we display the expected improvement C b /C s in complexity as a function of N for fixed F and K. We check for actual improvements with our implementation in the next section. The effect of the proposed method naturally increases linearly with the number of frequency bands considered.

EVALUATION

In this section, we show that the simplified structured covariance model [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF] leads to performance similar to the unconstrained PSDTF model [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF], while allowing for a significant computational speedup.

Experimental conditions

We consider the same evaluation material as presented in [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF], i.e., 4 synthetic mixtures, each of which being composed of K = 3 instrumental notes first presented separately, and then combined as various chords. Interestingly, the notes have a significant amount of overlap in the frequency domain. This dataset does not correspond to a real challenge in terms of modern source separation technology (see for instance the latest SiSEC report [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF] for that matter). However, they are sufficient for our purpose because our objective here is indeed only to assess whether the proposed structured model leads to a degradation in performance compared to PSDTF, as a price for its computational effectiveness. Comparing both methods on the same dataset seems sufficient for this.

In terms of metrics, we compute the Source to Distortion Ratio (SDR [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]) between the true source signals and the estimates, both for PSDTF and for the fast PSDTF we propose, abbreviated as fPSDTF here. We display the difference between these two scores for all mixtures as a function of the order N of the approximation in Fig. 4. Then, for one of the

Experimental results

As shown in these Figures, the fPSDTF has roughly the same performance than PSDTF, even for small N . As expected, the difference decreases with increasing N . This fact strongly supports the proposed model as a good approximation to PS-DTF that allows for effective implementations. Note that in our implementation, the computational speedup we observe was not as important as displayed on Fig. 3, but was rather of one order of magnitude, which is already noticeable. This is due to the fact that the PSDTF implementation benefited from multicore architectures, while fPSDTF did not.

Finally, for one particular excerpt and source k, we display the basis matrices |L k | learned with PSDTF for different values of N , along with the corresponding loading factors S k . We can see that only a few basis vectors turn out to be active in this example, suggesting that assuming a diagonal plus lowrank structure for spectral covariances fits music analysis.

CONCLUSION

In this paper, we introduced a low-rank plus diagonal structure for approximating the positive semidefinite covariance matrices that are routinely used in source separation studies. We showed how straightforward applications of classical linear algebra methods could allow the inversion of sums of such matrices, leading to considerable speedups as compared to naive inversion. Considering the particular positive semidefinite tensor factorization (PSDTF) model, where those covariances are rather large, we showed that this approximation does not lead to noticeable in performance as compared to the unconstrained model, while operating significantly faster. Interestingly, the proposed methodology and structured approximation may be used whenever large covariance matrices are considered, which may also happen in the case of massively multichannel signals. Future work includes using the potential of this structured approximation for scaling up the method to large-scale separation scenarios.

Fig.

  Fig. Diagonal plus low-rank approximation of a full covariance matrix over frequency bins. L k is F × N with N F . [v] is a diagonal matrix with diagonal v.

  P k ] and an independent deterministic component x (d) kt with low-rank covariance Y (d)

1 t

 1 kt . For this purpose, and for each component c = s or d, we need the respective Wiener gains W and the total posterior covariance as follows:for c ∈ {s, d} , Σ the re-estimation of h kt , we also need the total posterior covariance Σ kt of x kt , which is different from Σ kt are not independent conditionally on x t . It is straightforward to show that

Fig. 3 .

 3 Fig. 3. Theoretical speedup C b /C s over baseline PSDTF yield by the proposed inversion procedure permitted by the structured covariance model (10), for K = 5 sources and F = 1024 frequency bins. Observing this speedup in practice requires careful implementation of the inversion algorithm 2.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Performance increase yield by the proposed approximation fPSDTF over baseline PSDTF, for different values of the low-rank component order N .

Fig. 6 .

 6 Fig. 6. Basis matrices L k and the corresponding loading factors S k for: a) N = 5, b) N = 10 and c) N = 50. We can notice that even for large N , only a small number of basis are typically active.

Note that the computational complexity of the method proposed in[START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF] is comparable to that of the proposed EM algorithm. Even if it doesn't require computation of Σ kt , it requires Y -1 t V k , also scaling as O KT F 3 .

This work was partly supported by the research programmes KAMoulox (ANR-15-CE38-0003-01) funded by ANR, the French State agency for research.

† This work was partly supported ACCEL No. JPMJAC1602 and KAK-ENHI Nos. 26700020 and 16H01744 funded by JST, Japan.


B V Uspensky

K V Avramov

NONLINEAR NORMAL MODES OF STRONGLY NONLINEAR PERIODICALLY EXCITED PIECEWISE LINEAR SYSTEMS

A method for the numerical analysis of nonlinear normal modes of forced vibrations in strongly nonlinear systems with piecewise linear elastic characteristics is proposed. The approach is based on the combination of the Shaw-Pierre method of nonlinear normal modes with the Rauscher technique. As a result of application of this approach, the nonautonomous piecewise linear system is transformed into an autonomous system. For this system, we determine the Shaw-Pierre nonlinear normal modes. We also study the nonlinear torsional vibrations of the power transmission in a three-cylinder transport engine.

Introduction

The systems with finitely many degrees of freedom and piecewise linear elastic characteristics describe a broad class of engineering objects [1, pp. 59-74]. Numerous approaches were developed for the investigation of systems of this kind. Thus, in [START_REF] Wolf | Effect of smoothing piecewise-linear oscillators on their stability predictions[END_REF][START_REF] Begley | A comparison of piecewise linear and continuous approximating models[END_REF], it was shown that the approximation of piecewise linear characteristics by polynomial curves may lead to incorrect results in the study of stability of solutions and to significant errors in the numerical analysis of periodic motions of systems with strong nonlinearities. The analytic investigation of free oscillations in systems with piecewise linear elastic characteristic was carried out in [START_REF] Chati | Modal analysis of a cracked beam[END_REF][START_REF] Todd | Natural frequency considerations of an impact oscillator[END_REF][START_REF] Vestroni | A perturbation method for evaluating nonlinear normal modes of a piecewise linear twodegrees-of-freedom system[END_REF]. The procedure of finding Shaw-Pierre nonlinear normal modes of free oscillations in piecewise linear systems is presented in [START_REF] Jiang | Large-amplitude nonlinear normal modes of piecewise linear systems[END_REF]. The response of piecewise linear systems to the action of harmonic excitations was analyzed in [START_REF] Luo | The mapping dynamics of periodic motions for a three-piecewise linear system under a periodic excitation[END_REF][START_REF] Ji | Dynamics of a piecewise linear system subjected to a saturation constraint[END_REF][START_REF] Kim | Super-and sub-harmonic response calculations for a torsional system with clearance nonlinearity using the harmonic balance method[END_REF][START_REF] Zhu | Nonlinear dynamics of a one-way clutch in belt-pulley systems[END_REF]. The chaotic oscillations in piecewise linear systems with two and three degrees of freedom were studied in [START_REF] Avramov | Regions of chaotic oscillations of discrete mechanical systems with piecewise-linear elastic characteristics[END_REF][START_REF] Avramov | Resonance under random vibrations of discrete dynamic systems with piecewise-linear elastic characteristics[END_REF]. For the analysis of piecewise linear systems, it was proposed to use amplitude surfaces [START_REF] Avramov | Bifurcation analysis of a vibropercussion system by the method of amplitude surfaces[END_REF]. The nonlinear normal modes of free oscillations of piecewise linear systems were considered in [START_REF] Sh | Normal modes for piecewise linear vibratory systems[END_REF][START_REF] Jiang | Large-amplitude nonlinear normal modes of piecewise linear systems[END_REF].

There are only several works devoted to the numerical analysis of nonlinear normal modes under the conditions of forced vibrations. It seems likely that the Shaw-Pierre nonlinear normal modes under the conditions of forced vibrations were first studied in [START_REF] Jiang | Nonlinear normal modes for vibratory systems under harmonic excitation[END_REF]. The possibility of combination of the method of nonlinear normal modes with the Rauscher technique for the analysis of the Kauderer-Rosenberg nonlinear normal modes under the conditions of forced vibrations was studied in the book [19, pp. 230-238]. The forced and parametric vibrations were computed with the help of the method of nonlinear normal modes and the Rauscher technique in [START_REF] Avramov | Analysis of forced vibrations by nonlinear modes[END_REF][START_REF] Avramov | Nonlinear modes of parametric vibrations and their applications to beams dynamics[END_REF]. A detailed survey of works dealing with the theory of nonlinear normal modes can be found in [START_REF] Avramov | Review of applications of nonlinear normal modes for vibrating mechanical systems[END_REF][START_REF] Mikhlin | Nonlinear normal modes for vibrating mechanical systems. Review of theoretical developments[END_REF].

In the present work, forced vibrations in strongly nonlinear piecewise linear systems with any number of degrees of freedom are studied by the Shaw-Pierre method of nonlinear normal modes and the Rauscher technique. As a result of application of this approach, the nonautonomous dynamical systems are reduced to equivalent autonomous systems, which are used to compute the Shaw-Pierre nonlinear normal modes.

Statement of the Problem

We study a system with N degrees of freedom subjected to the action of external harmonic excitations. The motion of this system is analyzed in a system of generalized coordinates x i , i = 1,…, N :

!! x i + f i (x 1 ,…, x N ) = 0, i = 1,…, k -1, k + 1,…, N , !! x k + f k (x 1 ,…, x N ) = A cos ωt . (1)
The nonlinear functions f i , i = 1,…, N , are trilinear, depend on the vector of generalized coordinates x =

x 1 ,…, x N [ ] T , and take the form

f i (x) = k 1i x, Δ 13 < h T x < Δ 12 , k 2i x + b 2i , h T x ≥ Δ 12 , k 3i x + b 3i , h T x ≤ Δ 13 , ⎧ ⎨ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ (2)
where k 1i , k 2i , and k 3i are the row vectors of coefficients, b 2i and b 3i are free terms guaranteeing the con- tinuity of the elastic characteristic at the points of kinks, and h is the column vector. The aim of the present work is to determine the nonlinear normal modes under the conditions of forced vibrations in the strongly nonlinear system [START_REF] Avramov | Nonlinear Dynamics of Elastic Systems[END_REF]. The detailed description of the definitions and properties of nonlinear normal modes can be found in the book [START_REF] Avramov | Nonlinear Dynamics of Elastic Systems[END_REF] and in the works [START_REF] Pesheck | A new Galerkin-based approach for accurate nonlinear normal modes through invariant manifolds[END_REF][START_REF] Wolf | Effect of smoothing piecewise-linear oscillators on their stability predictions[END_REF][START_REF] Zhu | Nonlinear dynamics of a one-way clutch in belt-pulley systems[END_REF].

Numerical Analysis of Nonlinear Normal Modes for Nonautonomous Systems

We propose a procedure of finding forced vibrations in strongly nonlinear piecewise linear systems based on the Rauscher technique of reduction of nonautonomous dynamical systems to equivalent autonomous systems. In the obtained autonomous dynamical system, we seek the Shaw-Pierre nonlinear normal modes. In this way, we determine the nonlinear normal modes under the conditions of forced vibrations.

The dynamical system (1) without external action can be rewritten in the normal coordinates of one of its linear parts. As this part, we choose the first row of function [START_REF] Luo | The mapping dynamics of periodic motions for a three-piecewise linear system under a periodic excitation[END_REF]. Then the free oscillations in this part are described by the dynamical system

!! x + K 1 x = 0, (3)
where the matrix K 1 has the block structure:

K 1 = k 11 k 12 … k 1N [ ] T .
The blocks of this matrix are presented in Eq. ( 2). The dynamical system (3) has the vector of modal coordinates η = [η 1 ,…, η N ] connected with the vector x by the following nonsingular linear transformation of coordinates:

x = Qη,
where Q is the matrix of eigenvectors of the matrix K 1 . In the modal coordinates, the dynamical system (3) takes the form

!! η i + ω i 2 η i = 0 , i = 1,…, N ,
where ω i 2 are the eigenvalues of the matrix K 1 .

The dynamical system (1) can be represented in the modal coordinates η as follows:

!! η + Λ 1 η = " f (η) + Q * , k M k A cos ωt , (4)
where M k is the mass of the k th object in the system, Λ 1 is a diagonal matrix formed by the eigenvalues of the matrix

K 1 , Q * ,k = (Q 1,k ,…, Q N ,k
) is the k th column of the matrix Q , k is the number of mass subject- ed to the action of external periodic forces in system (1). The vector function ! f (η) can be represented in the form

! f (η) = 0, Δ 13 < h T Qη < Δ 12 , Ζ 2 η + ξ 2 , η ≥ Δ 12 , Ζ 3 η + ξ 3 , η ≤ Δ 13 , ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪
where

Ζ 2 = Λ 1 -Q -1 K 2 Q , Ζ 3 = Λ 1 -Q -1 K 3 Q , ξ 2 = Q -1 b 2 , ξ 3 = Q -1 b 3 , K 2 = k 21 k 22 … k 2 N ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ , K 3 = k 31 k 32 … k 3N ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ , b 2 = b 21 b 22 … b 2 N ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ , b 3 = b 31 b 32 … b 3N ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ .
First, we consider the nonlinear normal modes in the autonomous dynamical system that follows from (4) if the external periodic forces are neglected. According to the procedure of finding the Shaw-Pierre nonlinear normal modes [START_REF] Avramov | Nonlinear Dynamics of Elastic Systems[END_REF][START_REF] Avramov | Analysis of forced vibrations by nonlinear modes[END_REF], we choose one of the modal coordinates η and its velocity as driving. We denote these coordinates by (η 1 , ! η 1 ). All remaining coordinates of the vector η are regarded as driven. In this case, the driven coordinates of the nonlinear normal modes can be expressed via the driving coordinates in the following way:

η k = Θ k (η 1 , ! η 1 ) , ! η k = Ξ k (η 1 , ! η k = 2,…, N . (5)
We now perform a change of variables. As the driving coordinates, we choose coordinates a and ϕ con- nected with η 1 and ! η 1 by the formulas

η 1 = a cos ϕ , ! η 1 = -aω 1 sin ϕ. (6)
We introduce representation [START_REF] Avramov | Regions of chaotic oscillations of discrete mechanical systems with piecewise-linear elastic characteristics[END_REF] in the first relation in [START_REF] Avramov | Nonlinear modes of parametric vibrations and their applications to beams dynamics[END_REF]. The function obtained as a result can be represented in the form of truncated Fourier series as follows [START_REF] Avramov | Nonlinear Dynamics of Elastic Systems[END_REF]:

η i = P i (a, ϕ) = C i ℓ,m A ℓ (a) cos (m -1) ϕ m=1 Nϕ ∑ ℓ=1 N a ∑ , (7)
where A ℓ (a) are basis functions, N a is the number of basis functions A ℓ , and N ϕ is the number of elements of the expansion in ϕ.

We now study forced vibrations in the essentially nonlinear system (4). It is known [START_REF] Avramov | Nonlinear Dynamics of Elastic Systems[END_REF] that the forced vibrations in this system are well described in the harmonic approximation:

η i = γ i cos ωt . ( 8
)
Comparing relations ( 8) and ( 6), we get

A cos (ω t) = A cos ϕ = A η 1 η 1 2 + ! η 1 2 1 ω 1 2 . ( 9
)
Thus, the external force in the dynamical system (4) is expressed in terms of the generalized coordinate η 1 and velocity ! η 1 . Hence, the nonautonomous dynamical system (4) can be replaced by the so-called equivalent pseudoautonomous dynamical system. For this purpose, we insert relations ( 9) in ( 4) and get the following pseudoautonomous dynamical system:

!! η + Λ 1 η = g(η), ( 10
)
where the elements of the matrix g = [g 1 ,…, g N ] are given by the formula

g i (η) = f i (η) + AQ i,k η 1 M k η 1 2 + ! η 1 2 1 ω 1 2 .
According to [START_REF] Mikhlin | Nonlinear normal modes for vibrating mechanical systems. Review of theoretical developments[END_REF], the modal coordinates (a, ϕ) satisfy the system of ordinary differential equations:

! a = - g 1 (a, ϕ) ω 1 cos ϕ, ! ϕ = ω 1 - g 1 (a, ϕ) aω 1 sin ϕ. ( 11
)
The determination of the invariant manifolds of system ( 1) is reduced to finding the coefficients

C i ℓ,m , i = 2,…, N , ℓ = 1,…, N a , m = 1,…, N ϕ ,
of expansion [START_REF] Avramov | Resonance under random vibrations of discrete dynamic systems with piecewise-linear elastic characteristics[END_REF]. For this purpose, we represent the first and second derivatives in Eqs. [START_REF] Chati | Modal analysis of a cracked beam[END_REF] in the form

! η i = C i ℓ,m ′ A ℓ cos (m -1)ϕ ⋅ ! a -A ℓ (m -1) sin (m -1)ϕ ⋅ ! ϕ ( ) m=1 N ϕ ∑ ℓ=1 N a ∑ = C i ℓ,m g 1 (a, ϕ) ω 1 ⎛ ⎝ ⎜ -′ A ℓ cos ϕ cos (m -1) ϕ ( m=1 N ϕ ∑ ℓ=1 N a (m -1) A ℓ a sin ϕ sin (m -1)ϕ ⎞ ⎠ ⎟ -ω 1 (m -1)A ℓ sin (m -1)ϕ ⎞ ⎠ ⎟ , (12)
!! η i = C i ℓ,m ′′ A ℓ cos (m -1)ϕ ⋅ ! a 2 ( m=1 Nϕ ∑ ℓ=1 N a ∑ + ′ A ℓ cos (m -1)ϕ ⋅ !! a -2(m -1) ′ A ℓ sin (m -1)ϕ ⋅ ! a ! ϕ -A ℓ (m -1) 2 ! ϕ 2 cos (m -1)ϕ -A ℓ sin (m -1)ϕ ⋅ !! ϕ ) = C i ℓ,m cos (m -1)ϕ ′′ A ℓ ! a 2 + ′ A ℓ !! a ( ) ( m=1 Nϕ ∑ ℓ=1 N a ∑ -sin (m -1)ϕ 2(m -1) ′ A ℓ ⋅ ! a ! ϕ + A ℓ !! ϕ ( ) ) ,
where

′ A ℓ = dA ℓ da and ′′ A ℓ = d 2 A ℓ da 2 .
We now represent the second derivatives of modal coordinates in Eqs. [START_REF] Ji | Dynamics of a piecewise linear system subjected to a saturation constraint[END_REF] in the following form:

!! a = d dt - g 1 ω 1 sin ϕ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = - 1 ω 1 ∂g 1 ∂a ! a + ∂g 1 ∂ϕ ! ϕ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ sin ϕ + g 1 cos ϕ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ , ( 13
) !! ϕ = d dt ω 1 - g 1 aω 1 cos ϕ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = - 1 ω 1 ∂g 1 ∂a ! a + ∂g 1 ∂ϕ ! ϕ - g 1 a ! a ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 1 a cos ϕ - g 1 a ! ϕ sin ϕ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ , (14)
where

∂g 1 ∂a = 0, Δ 2 < h T Qη < Δ 3 , Ζ 2 1i ∂η i ∂a i=1 N ∑ , h T Qη ≤ Δ 2 , Ζ 3 1i ∂η i ∂a i=1 N ∑ , h T Qη ≥ Δ 3 , ⎧ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ∂g 1 ∂ϕ = - AQ i,k M k sin ϕ, Δ 2 ≤ h T Qη ≤ Δ Ζ 2 1i ∂η i ∂ϕ i=1 N ∑ - AQ i,k M k sin ϕ, h T Qη < Δ 2 , Ζ 3 1i ∂η i ∂ϕ i=1 N ∑ - AQ i,k M k sin ϕ, h T Qη > Δ 3 , ⎧ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ( 15
)
∂η i ∂a = cos ϕ, i = 1, C i ℓ,m ′ A ℓ cos (m -1)ϕ m=1 N ϕ ∑ ℓ=1 N a ∑ i ≠ 1, ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ ∂η i ∂ϕ = -a sin ϕ, i = 1, - C i ℓ,m (m -1)A ℓ sin (m -1)ϕ, i ≠ 1. m=1 N ϕ ∑ ℓ=1 N a ∑ ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ (16)
We now substitute relations ( 15) and ( 16) in ( 13) and ( 14) and then insert the results in Eq. [START_REF] Chati | Modal analysis of a cracked beam[END_REF]. Finally, we get the discrepancy vector for the solutions of the system of equations [START_REF] Chati | Modal analysis of a cracked beam[END_REF]:

Z(a, ϕ, C) = !! η(a, ϕ, C) + Λ 1 η(a, ϕ, C) -g(a, ϕ, C), ( 17
)
where C is a vector formed by the elements C i ℓ,m , i = 2,…, N , ℓ = 1,…, N a , m = 1,…, N ϕ . By virtue of relations ( 6) and ( 7), the first element of the discrepancy vector ( 17) is equal to zero. The remaining elements are nonzero. In order to find the nonlinear normal modes, we use the Galerkin method and write the conditions of orthogonality of the discrepancy vector Z(a, ϕ, C) to the basis functions A ℓ cos (m -1)ϕ:

(!! η(a, ϕ, C) + Λ 1 η(a, ϕ, C) -g(a, ϕ, C))A ℓ cos [(m -1)ϕ] da dϕ = 0 a 0 a 1 ∫ 0 2π ∫ . ( 18
)
Thus, relations (18) form a system of (N -1)N a N ϕ nonlinear algebraic equations for the elements of the vector C . They can be solved numerically by the Newton-Raphson method.

To study the on the invariant manifolds, we choose the initial values of the modal variables a(0) and ϕ(0). The motion on the form is determined by the numerical integration of the equations of modal dynam- ics [START_REF] Sh | Normal modes for piecewise linear vibratory systems[END_REF]. To analyze the stability of the obtained periodic motions, we numerically solve the equations in variations. We construct the corresponding fundamental matrix and determine its eigenvalues (multiplicators) [START_REF] Avramov | Nonlinear Dynamics of Elastic Systems[END_REF].

Numerical Analysis of Nonlinear Torsional Vibrations

We consider a model of torsional vibrations of the power transmissions in engines of internal combustion. The computational scheme of a system of this kind is presented in Fig. 1. We analyze the power transmission of a three-cylinder transport engine with oppositely moving pistons. The concentrated disks of the model describe the torsional vibrations of two crankshafts of the engine. The nonlinear elastic characteristic describes the final transmission between the crankshafts. It consists of five gears one of which contains an elastic clutch. The forced torsional vibrations of the power transmission are described by the following system of ordinary differential equations:

I 1 !! θ 1 + f (θ 1 -θ 2 ) = A cos (ωt), I 2 !! θ 2 -f (θ 1 -θ 2 ) + c 1 θ 2 = 0 , (19)
where

f (θ 1 -θ 2 ) = c 1 (θ 1 -θ 2 ), -Δ < θ 1 -θ 2 < Δ, c 2 (θ 1 -θ 2 ) + Δ(c 2 -c 1 ), θ 1 -θ 2 ≥ Δ, c 2 (θ 1 -θ 2 ) -Δ(c 2 -c 1 ), θ 1 -θ 2 ≤ -Δ, ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪
θ 1 and θ 2 are generalized coordinates used to describe the vibrations of two crankshafts of the engine, and c 1 is the stiffness of the gear transmission between the crankshafts. We performed the numerical analysis for the following dimensionless values of the parameters of the system:

c 1 = 300 , c 2 = 150 , Δ = 1.04 , I 1 = 0.3, I 2 = 0.45 , A = 77.7 .
The eigenfrequencies of linear low-amplitude vibrations of the system are as follows:

ω 1 = 18.25 and ω 2 = 44.72 .

The nonlinear normal modes under the conditions of forced vibrations were studied by using the approach described above. The numerically obtained nonlinear normal modes were rearranged in the configuration space (θ 1 , θ 2 ) . The results of numerical modeling of the nonlinear normal modes were studied for the ampli- tude of vibration a(0) = 20 Δ. The frequency of perturbing force was set equal to ω = 17.186 rad/sec. The re- sults of numerical analysis of this nonlinear normal mode performed by using the proposed algorithm are illustrated by the continuous line in Fig. 2. To confirm the obtained nonlinear normal mode, we carried out the direct numerical integration of the system of differential equations [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF] for 100 periods of action of the perturbing force. The results of calculations are presented in the form of squares in Fig. 2. It is easy to see that the obtained trajectories are close to each other, which confirms the validity of our calculations of the nonlinear normal mode.

We also computed the nonlinear normal modes in the region of the first main resonance of the system. The results are shown in Fig. 3, where the amplitude of vibrations θ 1 (t) is plotted on the ordinate and the frequency of perturbing force ω is plotted on the abscissa. The dash-dotted line reflects the skeleton curve of free vibrations, the solid line shows the branches of stable forced vibrations of the system, and the dashed line corresponds to the regions of unstable forced vibrations of the system. In order to verify the accumulated results, we performed the direct numerical integration of the dynamical system [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF] with initial conditions corresponding to the nonlinear normal mode of forced vibrations. The results of this direct integration are shown by symbols in the amplitude-frequency characteristic. It is easy to see that the data obtained by using both methods are quite close.

CONCLUSIONS

We propose a new approach to the determination of the Shaw-Pierre nonlinear normal modes of forced vibrations in strongly nonlinear systems with piecewise linear elastic characteristic. It is based on the combination of the method of Shaw-Pierre nonlinear normal modes for piecewise linear systems in the process of free vibration with the Rauscher technique.

We studied forced vibrations in a strongly nonlinear system with two degrees of freedom, which are nonlinear normal modes. They are presented in the amplitude-frequency characteristic, which has the ordinary form. Thus, the forced vibrations in strongly nonlinear piecewise linear systems can be represented in the form of nonlinear normal modes.

Fig. 1 .

 1 Fig. 1. Mechanical system.

Fig. 2 .

 2 Fig. 2. Nonlinear normal mode in the configuration space for ω = 17.186 rad/sec.

Fig. 3 .

 3 Fig. 3. Amplitude-frequency characteristic of the response.