domain: math nlin spi

text:
J.-P Giovanangeli

C Kharif

Y A Stepanyants

SOLITON SPECTRA OF RANDOM WATER WAVES IN SHALLOW BASINS

Keywords: Mathematics Subject Classification. 76B15, 76B25, 35Q51, 35Q53, 37K40 Shallow water, wind waves, random wave field, wave spectra, solitons, numerical modelling, laboratory experiments

Soliton spectra of random water waves in shallow basins

Introduction

The traditional approach in the problem of wind wave study is based on the analysis of Fourier spectra and determination of their peculiarities. There is a vast number of papers both theoretical and experimental where this problem has been considered; it is impossible to list all of them here. Therefore, we refer only to the review chapter "Wind waves" by Zaslavsky and Monin in the book [14] where a reader can find key references in this field. A lot of interesting and useful information has been obtained about wind waves, and their analysis and interpretation have been implemented in terms of Fourier spectra.

Meanwhile, the Fourier analysis of wind waves provides researchers with only some piece of objective information whereas many important features of wind waves remain hidden. One of the serious obstacles making the Fourier analysis ineffective in application to surface oceanic waves is the nonlinear character of such waves, whereas the Fourier analysis is a linear operation applicable to systems obeying the superposition principle.

Osborne with co-authors (see, e.g., [START_REF] Osborne | The inverse scattering transform: tools for the nonlinear Fourier analysis and filtering of ocean surface waves[END_REF][START_REF] Osborne | Solitons in the periodic Korteweg-de Vries equation, the Θ-function representation, and the analysis of nonlinear, stochastic wave trains[END_REF][START_REF] Osborne | Nonlinear Ocean Waves and the Inverse Scattering Transform[END_REF][START_REF] Osborne | Soliton basis states in shallow-water ocean surface waves[END_REF][START_REF] Segre | Analysis of laboratory surface wave data using the scattering transform for the periodic Korteweg-de Vries[END_REF] and references therein) have developed the method of nonlinear spectral analysis of shallow water waves described by the Korteweg-de Vries (KdV) equation. Osborne's approach is based on the application of the inverse scattering method (ISM) to the analysis of random field data in a one-dimensional space domain with the periodic boundary conditions. The main idea of his approach was in the presentation of complex initial disturbance in terms of a set of elliptic functions (cnoidal waves). These functions can be considered as the nonlinear eigenmodes which are preserved in the process of wave field evolution in contrast to the linear sinusoidal eigenmodes. This means that a nonlinear wave spectrum calculated on the basis of these modes is invariant in time while a usual Fourier spectrum is variable due to nonlinear interactions between the different sinusoidal harmonics. An important feature of nonlinear eigenmodes is their reducibility to the sinusoidal eigenmodes in the case of small amplitudes. In other words, a nonlinear spectrum naturally reduces to the Fourier spectrum if the analyzed wave field is quasi-linear. However, the mathematical and numerical machinery used for calculation of nonlinear eigenmodes is not simple in contrast to the linear case and, apparently, it is impractical especially for the applied mariner engineers.

Here we propose a very similar to Osborne's, but a bit different approach to the analysis of random water waves that is also based on the application of ISM. The essential feature of our approach is the interpretation of a random initial wave field in terms of an ensemble of solitons and quasi-linear ripples rather than the set of elliptic eigenmodes (a similar approach was recently realised in Ref. [START_REF] Costa | Soliton turbulence in shallow water ocean surface waves[END_REF]). The idea is illustrated by an example of shallow water waves described by the classical KdV equation with the random initial data. If one takes some portion of random field data (which should be long enough), the number of solitons, their amplitudes, speeds, characteristic durations, etc., can be calculated then by means of ISM (e.g., by solving numerically the associated Sturm-Louiville problem) or by the direct numerical simulation of the corresponding KdV equation. In both cases, the numerical codes are currently very well developed and easily available. In the meantime, the knowledge of number of solitons obscured in the random wave field, their parameters and statistics is a matter of independent interest per se. We describe our approach below in detail and give some examples (preliminary results were reported at the conference OCEANS'13 MTS/IEEE in San Diego, USA [START_REF] Giovanangeli | Nonlinear spectra of shallow water waves[END_REF]).

Before we start, it is useful to remind that the soliton turbulence of rarified soliton ensembles in strongly integrable systems is trivial to certain extent -the distribution function of solitons is unchanged in time [START_REF] Zakharov | Kinetic equation for solitons[END_REF][START_REF] Zakharov | Turbulence in integrable systems[END_REF] (the definition of strongly and weakly integrable systems is given in [START_REF] Zakharov | Turbulence in integrable systems[END_REF]). This is a consequence of trivial character of soliton interactions in such systems. The solitons do not change their parameters after collisions, and only the paired collisions occur between them. However, the dynamics of a dense ensemble of solitons is more complicated and soliton turbulence is nontrivial even in the integrable nonlinear wave equations [START_REF] El | The thermodynamic limit of the Whitham equations[END_REF][START_REF] El | Critical density of a soliton gas[END_REF][START_REF] El | Kinetic equation for a dense soliton gas[END_REF]. In particular, in reference [START_REF] El | Critical density of a soliton gas[END_REF] the quantitative criterion of the term "dense soliton gas" was introduced and was shown that the density of KdV soliton gas is bounded from above. The critical gas density, apparently, depends on soliton distribution function, which makes important the determination of such function in a concrete physical problem such as water wave turbulence. Numerical experiments confirming the developed theory in reference [START_REF] El | Critical density of a soliton gas[END_REF] for the particular model distribution functions were reported in [START_REF] Carbone | Macroscopic dynamics of incoherent soliton ensembles: soliton gas kinetics and direct numerical modelling[END_REF].

The KdV model considered here is the typical example of strongly integrable system applicable to real physical systems. This makes topical the development of handy methods of extraction of soliton distribution function from natural complex fields. One of the experimental approaches to the solution of this practical problem has been considered for internal waves in Okhotsk Sea [START_REF] Nagovitsyn | Observation and analysis of solitary internal waves in the coastal zone of the Sea of Okhotsk. (English translation of the Russian journal Morskoi Gidrofizicheskii Zhurnal)[END_REF] and another example of processing of surface wave observational data was reported in reference [START_REF] Costa | Soliton turbulence in shallow water ocean surface waves[END_REF]. We hope, this publication will stimulate further interest to this important problem.

The Korteweg-de Vries model and data processing

Let us assume that there is a data of recorded surface waves at some fixed point x 0 of a shallow-water basin. So that the elevation η of the water level at this point is the known function of time: η(x 0 , t) = f (t) where f (t) is some random function. A typical example of random surface waves usually measured in shallow-water basin is shown in Figure 1. Here these data were artificially produced by means of a computer using the random number generator just to illustrate the idea of our approach. For the description of further space evolution of data measured at the point x 0 , the so-called, timelike KdV (TKdV) equation [START_REF] Osborne | Soliton basis states in shallow-water ocean surface waves[END_REF] is used:

∂η ∂x + 1 c 0 ∂η ∂t -αη ∂η ∂t -β ∂ 3 η ∂t 3 = 0, (2.1)
where c 0 = gh; α = 3/(2c 0 h); β = h 2 /(6c 3 0 ) with g being the acceleration due to gravity and h being an unperturbed water depth.

In the process of evolution of an initial perturbation one can expect emergence of a number of solitons with different amplitudes and phases (e.g., time markers of soliton maxima in the wave record at a given point of observation). Soliton solution to the TKdV equation (2.1) has the form:

η(x, t) = A sech 2 t -x/V T , (2.2)
where the velocity V and duration T are relate to the amplitude A:

V = c 0 1 -c 0 αA/3 ≈ c 0 1 + c 0 αA 3 , T = 12β αA , (2.3)
the approximate formula is valid for small amplitude solitons when c 0 αA/3 1. As it was mentioned in [START_REF] Osborne | Soliton basis states in shallow-water ocean surface waves[END_REF], "it may not be possible to observe solitons in real space" because of numerous nonlinear interactions of solitons both with each other and with chaotic radiation background component. However, "it would be naive to conclude that solitons are not present or that their dynamics are not important" in the evolution of the initial perturbation. Inasmuch as solitons are very stable with respect to interaction with others wave perturbations and influence of external effects (such as viscosity, inhomogeneity, etc.), it is a matter of interest to extract them from the irregular components of a wave field and to describe their statistical properties and contribution to the total wave energy.

The solution of this problem can be done by the following way. Let us consider a very long portion of recorded measurement data of surface perturbation at any given point x 0 . The characteristic duration of this portion T p is assumed to be much greater than the typical soliton time scale T . Let us represent a perturbation with the help of some dimensionless function ϕ(t): η(0, t) = U ϕ(t/T p ), (2.4) where U is the characteristic wave "amplitude", e.g., the maximum value of perturbation η(0, t) in the considered portion of data. By means of the transformation

u = η U , ξ = - αU T p x, τ = 1 T p t - x c 0 . (2.5)
Equation (2.1) and the corresponding "initial" perturbation (2.4) can be reduced to the standard form [START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF] (the term "initial" used here as it is traditionally used in mathematics for the solution of Cauchy problem of differential equations, but in fact, the perturbation is given at the fixed spatial point, i.e., it is rather the boundary condition):

u ξ + uu τ + 1 σ 2 u τ τ τ = 0 (2.6) u(0, τ ) = ϕ(τ ) (2.7)
with one dimensionless parameter σ 2 known in the oceanography as the Ursell parameter and defined as

σ 2 = αU T 2 p β . (2.8)
As it was mentioned above, we consider the case when the duration of the perturbation is long enough, so that σ 2 1. In this case the number of solitons obscured in the "initial" perturbation is also very big in general, and it is reasonable to describe them by the distribution function f (A). This function determines the number of solitons dN within the interval (A, A + dA) [START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF] dN = f (A)dA.

(2.9)

According to the theory developed in [START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF], the distribution function can be calculated at large values of σ by means of the formula:

f (A) = σ 4π √ 3U L dτ 2U ϕ(τ ) -A , (2.10)
where the interval of integration L is determined by the condition 2U ϕ(τ ) > A.

(2.11)

As follows from equation (2.10), soliton amplitudes are distributed in the interval [START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF] 0

< A < 2U max[ϕ(τ )], (2.12)
and their total number can be found from the formula

N = ∞ 0 f (A) dA = σ π √ 6 ϕ(τ )>0 ϕ(τ ) dτ. (2.13)
Therefore for large σ the total number of solitons is determined only by those intervals of τ -axis where function ϕ(τ ) is nonnegative! As well known, the KdV equation possesses an infinite number of conserved densities I n [START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF][START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF][START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. One of them,

I 2 = +∞ -∞ η 2 (x, t) dt, ( 2.14)
is proportional to wave energy and therefore is of a special physical interest. For the sake of simplicity we will call this value simply the "energy". The fraction of energy of a nonsoliton component of a perturbation to the total energy of "initial" perturbation can be determined by means of formulae (18.27), (18.28) and (19.8) in reference [START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF]:

I ns 2 I tot 2 = η(0,τ )<0 η 2 (0, τ ) dτ Tp 0 η 2 (0, τ ) dτ . ( 2.15)
The total energy of a soliton component in the wave field can be readily calculated, if soliton amplitudes are known:

I sol 2 = 4 β 3α N k=1 A 3/2 k = 4h 3 √ 3g N k=1 A 3/2 k .
(2.16)

In the last expression the values of coefficients α and β for surface water waves were used (see above after Eq. (2.1)).

To illustrate the idea of this approach we consider below several examples of different "initial" perturbations.

Sinusoidal perturbation and its modifications

A sinusoidal perturbation

Let us study the KdV equation (2.1) with the periodical conditions in time. Let us set h = 1 m, then c 0 = 3.13 m/s, α = 0.479 s/m 2 , β = 0.00543 s 3 /m. Assume then that the perturbation has the form

η(0, t) = A sin ωt, (2.17)
where A = 0.01 m, ω = 0.0628 1/s. The value of Ursell parameter for this perturbation is σ 2 = 555 σ 2 s , where σ 2 s ≡ 12 is the value of Ursell parameter for a single soliton regardless of its amplitude and duration [START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF].

Figure 2 shows this sinusoidal initial perturbation and the result of its evolution at the distance x = 1.1 • 10 3 m when it is completely disintegrated into the sequence of solitons. After that moment, solitons begin interact with each other and the wave field looks much more complex. In Figure 2 one can distinguish nine different Table 1. Amplitudes of solitons, A k , emerging from the different initial perturbations. The fifth raw in the Table shows the relative difference RD in the soliton amplitudes of rows three and four in percents, this will be explained in Section 2.1.3.

Initial pert.

A 1. Soliton amplitudes were conditionally defined here as the difference between their maxima and the mean value of two nearest minima. The total energy of the sinusoidal initial perturbation over a period T = 2π/ω can be easily evaluated

k A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 8 A 9 (mm) (mm) (mm) (mm) (mm) (mm) (mm) (mm) (
I tot 2 = 1 2 T 0 η 2 (0, t) dt = A 2 2 T 0 sin 2 ωt dt = A 2 T 4 = 2.5 • 10 -3 m 2 s. (2.18)
Similarly, the energy of the non-soliton component of the perturbation is (see Eq. (2.15)) Thus, I ns 2 = I tot 2 /2, and the energy of all solitons in this case should be I sol 2 = I tot 2 /2 = 1.25 • 10 -3 m 2 s. Let us calculate however this energy directly by means of equation (2.16) using data of Table 1 for the soliton amplitudes

I ns 2 = η(0,τ )<0 η 2 (0, τ ) dτ = A 2 2 T T /2 sin 2 ωt dt = A 2 T 8 = 1.25 • 10 -3 m 2 s. (2.19)
I sol 2 = 4h 3 √ 3g 9 k=1 A 3/2 k = 4.922 • 10 -3 m 2 s. (2.20)
This value about four times greater then the expected. The discrepancy can be explained by the incorrect counting of number of solitons, as well as their amplitudes. As has been shown in [START_REF] Salupere | Korteweg-de Vries soliton detection from a harmonic input[END_REF][START_REF] Salupere | Long-time behaviour of soliton ensembles, part I -emergence of ensembles[END_REF][START_REF] Salupere | Long-time behaviour of soliton ensembles, part II -periodical patterns of trajectories[END_REF], the problem of detection of solitons emerging from the harmonic perturbations is not so simple even within the KdV equation (disintegration of sine wave onto set of solitons was recently studied also within the Gardner equation [START_REF] Kurkina | Nonlinear disintegration of sine wave in the framework of the Gardner equation[END_REF][START_REF] Kurkina | Long sine wave transformation in the framework of Gardner equation: spectral and statistical analysis[END_REF]). Actual number of solitons is always greater then that at the instant of time when they appear from the initial perturbation for the first time in the ordered form. Their amplitudes are also different from those which are seen in Figure 2. We will come back to this issue a bit later, and now let us consider a pulse-type initial perturbations defined on a compact support.

A pulse of sinusoidal profile

Consider now a pulse-type initial perturbation having the shape of one period of sine with the same amplitude and characteristic duration as in the previous case (see dashed line in Fig. 3). In the process of evolution this perturbation disintegrates into a sequence of five solitons (see solid line in Fig. 3), whose amplitudes are indicated in the third row of Table 1. An intense oscillatory tail behind the solitons is also appeared; the head portion of this tail is shown in Figure 3.

By comparison of rows two and three of Table 1, one can see that the number of solitons and their amplitudes are absolutely different in the periodic and nonperiodic cases. As all these solitons are practically independent at the distance indicated in Figure 3, their energies can be calculated independently. The calculation of cumulative energy of these five solitons yield:

I sol 2 = 4h 3 √ 3g 5 k=1 A 3/2 k = 1.2531 • 10 -3 m 2 s. (2.21)
This result already agrees quite well with the theoretical prediction.

A half-sine pulse

Let us consider now another pulse-type perturbation which contains only a half period of a sine-function of a positive polarity (see dashed line in Fig. 4).

This perturbation disintegrates into the same number of solitons as in Section 2.1.2. At the same distance x = 6 • 10 4 m, the soliton amplitudes are practically the same as in the previous case (cf. rows three and four in Tab. 1). The fifth raw in Table 1 shows the relative difference in the soliton amplitudes in percents, RD = (A i -A i )/A i , where A i and A i are the amplitudes of i-th solitons emerged from the sinusoidal and half-sine pulses respectively.

As one can see from Table 1 and Figures 3 and4, there is no much differences in the corresponding soliton amplitudes for the sinusoidal and half-sine pulses. In the next Subsection we will consider a simple analytical example to justify this numerical observation.

Analytical solutions to the associated Schrödinger equation for the rectangular and meander-type pulses

To understand better the regularity of distribution of amplitudes of solitons emerging from initial perturbations, let us consider two model initial perturbations: the rectangular pulse of positive polarity and amplitude U 0 (Fig. 5a) and meander-type pulse whose positive part coincides with the above rectangular pulse and negative part is the same but of opposite polarity (Fig. 5b). Assume that the duration of positive rectangular pulses are T .

According to the ISM [START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF][START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF][START_REF] Whitham | Linear and Nonlinear Waves[END_REF], to study the evolution of the initial perturbation within the KdV equation (2.1), one have to solve the complementary Schrödinger equation:

ψ (ξ) + σ 2 6 [ϕ(ξ) + E n ] ψ(ξ) = 0, (2.22)
where ψ is the auxiliary function; ϕ is the dimensionless function describing the shape of the initial perturbation (see above), it has the unit amplitude, i.e., ϕ max = 1, and the unit characteristic duration in the dimensionless variables (2.7). For the case a) in Figure 5, ϕ = 1 at 0 ≤ ξ ≤ 1 and zero beyond this interval. Similarly, for the case b) in Figure 5, ϕ = 1 at 0 ≤ ξ ≤ 1, ϕ = -1 at 1 < ξ ≤ 2 and zero beyond these intervals. Further, σ 2 is the Ursell parameter as defined in equation (2.10) with U ≡ U 0 and T p ≡ T ; E n < 0 with n = 1, 2, . . . , N are eigenvalues of the Schrödinger equation (2.22). The amplitudes of emerging solitons, A n , are related with the eigenvalues E n by the simple relation:

A n = -2U 0 E n . (2.23)
The analytical solution to the Schrödinger equation (2.22) with the rectangular potential function can be readily constructed (see e.g., [START_REF] Flügge | Practical Quantum Mechanics[END_REF][START_REF] Landau | Quantum Mechanics -Nonrelativistic Theory[END_REF]). Omitting simple, but tedious calculations, the result can be presented in the form of two transcendental equations determining eigenvalues E n of the Schrödinger equation (2.22):

tan σ 2 β 24α (1 + E) = -E 1 + E , (2.24) tan σ 2 β 24α (1 + E) = - 1 + E -E . (2.25)
In a similar manner the analytical solution can be found for the initial perturbation shown in Figure 5b. The outcome is presented by the following transcendental equation: where Th = tanh σ 2 β 6α (1 -E) .

tan σ 2 β 6α (1 + E) = - √ 1 -E 2 √ -E + √ 1 + E + √ 1 + E -E(1 + E) + 1 -E Th √ 1 -E 2 √ -E - √ 1 + E + -E(1 -E) - √ 1 + E(1 + E) Th , ( 2
Solutions of the transcendental equations (2.24)-(2.26) can be presented graphically as shown in Figures 6a andb, where the series of tan-type curves represent the tangential functions in the left-hand side of these equations, while line 1 represents the right-hand side of equation ( 2 2.

The last row in the table shows the relative difference in percents between the corresponding soliton amplitudes for the rectangular A r k and meander-type

A m k perturbations: RD = (A r k -A m k )/A r k .
As one can see from this table, the difference between the corresponding amplitudes A k is fairly small especially for the first largest solitons, and only for two last solitons of very small amplitudes the difference amounts for about 7% and 31% correspondingly. It can be readily shown that the larger the Ursell parameter, the greater the number of emerging solitons and the smaller the difference in theirs amplitudes.

Thus, one can conclude that the asymptotic theory developed in [START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF] can provide a good basis for the calculation of statistical properties of solitons obscured in the random wave field. Apparently, the most energetic Table 2. Amplitudes of solitons emerging from the rectangular and meander-type initial perturbations shown in Figure 5.

Initial pert.

A part of the soliton spectrum (a right wing of soliton distribution function at large amplitudes) is described fairly good while at small amplitudes the distribution function may be not quite correct.

Cosine initial pulse

Let us now consider why the number of solitons emerging in each period of a pure periodic perturbation (as well as their amplitudes) differs from those emerging from the pulse-type perturbations of the same shape, amplitude and duration (cf. Fig. 2 with Figs. 3 and4)? The answer is, apparently, as follows: In the pure periodic case the zero level of the physical system is not determined clearly. The system admits the minimum of the wave field as the zero level. Hence, a periodic sinusoidal perturbation can be treated as a periodic sequence of positive cosine-type pulses with respect to the minimum possible level. In a process of evolution, each pulse disintegrates into a number of solitons whose total amount is much greater than for sine-type perturbations considered above. New numerical simulation was carried out with this modified zero level for single cosine-type pulse as the initial perturbation; the result is shown in Figure 7.

There were clearly detected 11 solitons whose amplitudes are presented in Table 3. The amplitudes of solitons emerging from the periodic sinusoidal perturbation described above (see Sect. 2.1.1) are also shown in the same table for the comparison. This time soliton amplitudes for the sinusoidal perturbation were calculated more thoroughly using trace method suggested in [START_REF] Salupere | Korteweg-de Vries soliton detection from a harmonic input[END_REF]. When the sinusoidal perturbation disintegrates for the first time into the sequence of solitons ordered by their amplitudes as shown in Figure 2, the observer cannot see all Table 3. Amplitudes of solitons, A k , emerging from the cosine initial pulse and from the periodic sinusoidal perturbation.

Initial pert.

A

k A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 8 A 9 A 10 A 11 (mm) (mm) (mm) (mm) (mm) (mm) (mm) (mm) (mm) (mm) (mm)
Cos., pulse 37.75 33. [START_REF] Salupere | Long-time behaviour of soliton ensembles, part I -emergence of ensembles[END_REF] solitons because some of them are still obscured. We followed up for the development of the wave field after the distance x = 1.1 • 10 4 m which corresponds to Figure 2 and discovered that in the process of soliton interactions some more solitons of small amplitudes appear. So that the total number of solitons in this case was exactly the same as for the cosine initial pulse, N = 11. Their amplitudes were measured with respect to the minima of the initial sinusoidal perturbation. As one can see, the relative difference in percents between the corresponding soliton amplitudes for the cosine initial pulse and periodic sinusoidal perturbation, RD = (A c k -A s k )/A c k is not too big now and not exceeds 17%. (see the last row of Tab. 3).

Thus, if the zero level of the considered physical system is known and the perturbation eventually vanishes at the infinity then the statistics of obscured solitons is determined by positive humps of the perturbation with respect to this zero level. Namely such a situation takes place in the case of hydrophysical measurements in laboratory (water tanks) or in natural marine conditions. But if the perturbation is periodic in principle, then solitons, their numbers and amplitudes are determined by each hump of the perturbation with respect to the total minimum of the initial perturbation.

Data processing: a model example

A model spectrum and the range of its validity

Let us apply now the developed approach to the random perturbation artificially generated and presented in Figure 1. This perturbation was obtained by means of inverse Fourier transform of a series of 128 harmonics having random phases in the interval [0, 2π] and amplitudes distributed in accordance with the formula

S(ω) = S 0 (ω + ω 0 ) -2 tanh ω δω , (2.27)
where

S 0 = 2 • 10 -4 m, ω 0 = 2π • 10 -2 s -1 , and δω = 2.5 • 10 -2 s -1 .
For the construction of quasi-random wave field it was taken only the low-frequency and most energetic part of this spectrum, 0 ≤ ω ≤ ω lim , where ω lim = 0.767 s -1 (see vertical dashed line in Fig. 8). The limiting frequency, ω lim , was chosen for the following reasons. The dispersion relation for infinitesimal perturbations within the TKdV equation (2.1) is

k = ω c 0 + h 2 6c 3 0 ω 3 = ω c 0 1 + h 2 6c 2 0 ω 2 , (2.28)
where ω is the frequency of a sinusoidal wave, and k is the wave number, η ∼ e i(ωt-kx) . As is well known, this dispersion relation represents an approximation of a real dispersion relation of a physical system, e.g., water or plasma waves, when ω, k → 0. The range of validity of the dispersion relation (2.28) is restricted by the requirement that the second term in brackets is small in comparison with one (see, e.g., [START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF][START_REF] Karpman | Nonlinear Waves in Dispersive Media[END_REF][START_REF] Whitham | Linear and Nonlinear Waves[END_REF]). This condition gives

ω ω cr ≡ c 0 h √ 6 = 6g h . ( 2

.29)

In our case ω cr ≈ 7.67 s -1 . To satisfy condition (2.29), we set ω lim = 0.1ω cr = 0.767 s -1 .

Data processing of the model initial perturbation

Let us consider now a quasi-random wave field shown in Figure 1 and apply the approach developed above. Each positive hump enumerated in the figure can be studied separately by means of the numerical code for the TKdV equation (2.1). The numerical code was based on the explicit finite-difference scheme of a second-order accuracy both on spatial and temporal variables [START_REF] Yu | Modelling Nonlinear Wave Processes[END_REF]. The theoretical analysis shows that used central-difference scheme is conditionally stable provided that ∆x ∼ ∆t 3 , where ∆x and ∆t are the spatial and temporal mesh steps. The code works very effectively and fast so that the result of pulse fission on solitons was obtained very quickly. A typical picture of disintegration of one of the pulses (pulse No. 5 in Fig. 1) is shown in Figure 9.

The number of emerged solitons and their parameters can be easily calculated. To avoid errors in determination of soliton parameters, the numerical calculations were conducted until each soliton was sufficiently separated from others, so that their fields were not overlapped in the vicinity of their maxima. This procedure was carried out for each pulse shown in Figure 1. As the result, it was obtained a large number of solitons of different amplitudes, their total number for all pulses was N s = 73 which is enough for the illustrative purposes. The solitons were collected into several groups (15 groups in total) according to their amplitudes. This allowed us to construct a histogram of numbers of solitons versus their amplitudes. The histogram can be considered as a model of the distribution function f (A) of density of number of solitons on amplitudes (see Eq. (2.9)). The histogram obtained and the model distribution function for the considered illustrative example are shown in Figure 10. Frame (b) in the figure shows the cumulative distribution function, i.e., the total number of solitons whose amplitudes are not greater than the given value.

Laboratory experiment with wind waves on shallow water

The theory developed above was applied to the data of laboratory experiments with wind wave generation. Series of experiments were carried out in the Luminy (Marseilles) small tank made of plexiglass and having the following sizes (length × width × height): 865 cm × 64 cm × 50 cm. The water depth in the tank in different experiments ranged from 1 cm to 8 cm: h = 1, 2, 3, 4, 6 and 8 cm. Surface waves were generated by a wind blowing up over the water surface with the different mean velocities: V w = 5.29, 6.45, 8.62 and 13.24 m/s. The ventilator producing an air flow was installed at one of the ends of the tank. The blower width was the same as the width of the tank, 64 cm, but its height was 31 cm above the water level. The tank was covered above the blower by a plexiglass lid. At the opposite end of the tank it was placed a wave absorber to exclude reflected waves (an inclined bottom causing surface wave breaking). Two sensitive electric probes recording water level were placed at the distances 100 cm and 300 cm from the ventilator. The probes were thoroughly calibrated before each experiment; their sensitivities were equal and amounted s = 0.61 V/cm.

Below we present the analysis of only one of the series of experiments with the water depth h = 1 cm and wind velocity V w = 5.29 m/s. Other experimental series were analysed in a similar way. We have to make a reservation in advance that the experimental data with wind waves are not perfect for the illustration of suggested approach. Wind waves generated by permanently blowing wind is an active system, i.e., the system with the permanent energy pumping at each point of water surface. Moreover, a distributed external force due to wind is not a constant, it varies from some maximum value near the ventilator to some small value at the opposite end of the tank. This results in the different soliton distribution functions measured at two different distances from the ventilator. A small water viscosity can also affect the soliton distribution function.

Another difficulty with surface waves generated by wind is the effective generation of high frequency Fourier components, so that the essential portion of wave energy is contained in that part of Fourier spectrum which is beyond the range of validity of KdV equation. Therefore we were forced to restrict our analysis by only the lowfrequency components of the wave spectrum. Figure 11 presents the Fourier spectrum of wind waves recorded at two distances from the ventilator as indicated above. The critical frequency at h = 1 cm is ω cr = 76.7 s -1 , and to satisfy the condition (2.29) we cut the Fourier spectra of wind waves at ω lim = 13.4 s -1 (see dashed vertical line in Fig. 11) and ignored the high frequency portions of spectra above ω lim . The surface perturbation was reconstructed by means of the inverse Fourier transform on the basis of extracted portions of wave spectra in the range 0 < ω ≤ ω lim ). The fragments of 60-second duration records corresponding to two spatial points of measurements are shown in Figure 12.

As it is clearly seen from the comparison of two data records presented in frames (a) and (b), there is the essential difference between them, especially in the intensity of wave fields. This can be explained by the effect of a fetch on wind wave generation. Therefore, one can expect that the statistics of obscured solitons in the record of frame (b) is essentially richer than in frame (a).

The recorded data shown in Figure 12 were used as the input data for the TKdV equation (2.1). In the statistically equilibrium state each 60 s portion of recorded data is equivalent to the same portion taken at a different time, therefore one can expect that the number of solitons obscured in each portion of data is the same in average and their distribution function is invariant with respect to time shift. This was confirmed in the data processing.

The TKdV equation (2.1) was solved numerically using the recorded data of 60 s duration from the total time interval of 208 s. After a while solitons emerged from the quasi-random data, and their amplitudes were easily determined with the help of a special subroutine. This allowed us to determine the histogram of soliton numbers in the each particular interval of amplitudes A + ∆A; this is the analog of a differential distribution function. On the basis of this function we determined also the integral (cumulative) distribution function -the total number of solitons with the amplitudes less than A normalised by the total number of all solitons. Figure 13 The corresponding integral distribution functions for the experimental data of Figure 12 are shown in Figure 14 (lines 1) together with the approximative Poisson integral functions with the same parameters as in Figure 13. The total number of solitons emerged from the wave field shown in Figure 12a was 60, and emerged from the wave field shown in Figure 12b was 86. As expected, the time series recorded closer to the ventilator (Fig. 12a) contained less number of solitons than the time series recorded further from the ventilator (Fig. 12b). In the latter case the wave field was much better developed due to the influence of wave fetch.

If we assume that all 60 solitons in the time series shown if Figure 13a are uniformly distributed in the time interval of 60 s, then we obtain that the time interval per each soliton is ∆T 1 = 1 s. As follows from the histogram shown in Figure 13a the maximal number of solitons have amplitudes A m1 = 0.02 cm and the duration T s1 = 0.026 s. Therefore ∆T 1 /T s1 ≈ 3.85. The similar estimates for the time series shown if Figure 13b give the time interval per each soliton ∆T 2 ≈ 0.7 s. As follows from the histogram shown in Figure 13b the maximal number of solitons in this time series have amplitudes A m2 = 0.05 cm and the duration T s2 = 0.165 s. Therefore ∆T 2 /T s2 ≈ 4.2. Thus, we see that in both cases the "soliton gas" is very dense (cf. [START_REF] El | Critical density of a soliton gas[END_REF][START_REF] Shurgalina | Dynamics of irregular wave ensembles in the coastal zone[END_REF]). The fragments of numerical calculations presented in Figure 15 illustrate the soliton gas density in both time series.

Conclusion

To analyze long random time series of water waves in shallow basins we have proposed an approach which differs from the traditionally used Fourier analysis. Our approach is based on the extraction of obscured solitons from the complex wave fields and construction of histograms of solitons at different points of observation. The histograms can be considered as experimental counterpart of distribution functions of number of solitons on their amplitudes. According to the theoretical conception, a soliton component of a wave field in the well-developed nonlinear perturbations should dominate. As is well known, the number and individual parameters of solitons preserve in the conservative statistically homogeneous systems [START_REF] Zakharov | Kinetic equation for solitons[END_REF][START_REF] Zakharov | Turbulence in integrable systems[END_REF], therefore the distribution function (or histograms of solitons) is the same at different points of observation if the dissipative factors (i.e., viscosity or external sources of energy) are negligible. In contrast to that the Fourier spectrum changes due to nonlinearity.

Our approach is in line with the contemporary development of the theory of strong turbulence in the integrable or near-integrable systems [3, 5-9, 28, 31]. Experimentally constructed distribution function can be used for the determination of degree of soliton gas density -how far the density is from the critical value as defined in references [START_REF] El | Critical density of a soliton gas[END_REF][START_REF] Shurgalina | Dynamics of irregular wave ensembles in the coastal zone[END_REF]. Data processing of laboratory experiments presented in our paper supplement the data processing of field experiments reported in reference [START_REF] Costa | Soliton turbulence in shallow water ocean surface waves[END_REF].

A small dissipation can cause a gradual decay of soliton histograms and their distortion, in general. The histogram decay and its distortion depend on the specific type of dissipation; this problem has not been studied yet, although the decay of individual solitons under the influence of various types of dissipation has been investigated for the KdV [START_REF] Grimshaw | Damping of large-amplitude solitary waves[END_REF] and Benjamin-Ono [START_REF] Grimshaw | Decay of Benjamin-Ono solitons under the influence of dissipation[END_REF] solitons, as well as for the Kadomtsev-Petviashvili lumps [START_REF] Clarke | Decay of Kadomtsev-Petviashvili lumps in dissipative media[END_REF].

Our approach can provide some additional valuable information about the energy distribution in natural wave fields such as the relationship between the soliton and nonsoliton components of the perturbation, and may indicate on the existence and intensity of external sources or sinks of energy.

To determine the soliton number and amplitudes from the random time series, we applied direct numerical modelling for the evolution of initial data within the framework of the TKdV equation (2.1). The existing numerical codes (see, e.g., [START_REF] Yu | Modelling Nonlinear Wave Processes[END_REF]) allow us to obtain the results fairly quickly in the form convenient to further analysis. However, it is not the only method which can be applied; the numerical solution of the eigenvalue problem (2.22) can also be convenient and useful. Our experience with the application of the approach developed here shows that there is no problem with the determination of number and amplitudes of intense solitons, i.e., solitons of big and moderate amplitudes. However, it takes more efforts to determine parameters of solitons whose amplitudes are very small. In general, an uncertainty in the determination of parameters of solitons of very small amplitudes is higher than of moderate and big solitons. Meanwhile, the model example of Section 2.4.2 shown in Figure 10, as well as the experimental laboratory data shown in Figure 13 demonstarte that the number of such small-amplitude solitons may be relatively big.

As has been mentioned, the soliton distribution function remains unchanged in the integrable systems. However, the wave field randomly fluctuates in the process of evolution. This leads to the random fluctuations of local wave extrema. As has been shown in reference [START_REF] Shurgalina | Dynamics of irregular wave ensembles in the coastal zone[END_REF], the distribution function of wave extrema varies with time even when the wave field consists of solitons only. This can be of interest from the viewpoint of physical applications, but beyond the scope of this paper.

In conclusion, we emphasize that the approach developed here is applicable to the KdV-type systems, e.g., shallow-water waves (see, for example, Ref. [START_REF] Dytukh | Numerical simulation of a solitonic gas in KdV and KdV-BBM equations[END_REF] where the turbulence of soliton gas was studied both within the integrable KdV and non-integrable KdV-BBM equations). Its generalization to deep-water waves described by the Benjamin-Ono or nonlinear Shrödinger equation is the interesting and challenging problem.

Figure 1 .

 1 Figure 1. A typical example of random surface waves in a shallow water basin.

Figure 2 .

 2 Figure 2. Sinusoidal initial perturbation (dashed line) and the result of its evolution at the distance x = 1.1 • 10 4 m (solid line).

Figure 3 .

 3 Figure 3. A sinusoidal-shape pulse as the initial perturbation (dashed line) and the result of its evolution at the distance x = 6 • 10 4 m (solid line). Only a fragment of a very long quasi-sinusoidal tail behind five solitons is shown in the figure.

Figure 4 .

 4 Figure 4. A half-sine pulse as the initial perturbation (dashed line) and the result of its evolution at the distance x = 6 • 10 4 m (solid line).

Figure 5 .

 5 Figure 5. Two model initial perturbations: the rectangular pulse (a) and meander-type pulse with the zero mean value (b).

. 26 )Figure 6 .

 266 Figure 6. (a) Graphical solutions of transcendental equations (2.24) and (2.25). (b) Graphical solutions of transcendental equation (2.26).

  .24), line 2 represents the right-hand side of equation (2.25), and two branches of line 3 represent the right-hand side of equations (2.26). The plots were generated for the same values of Ursell parameter σ 2 = 555, the water depth h = 1 m, and the amplitude of initial perturbations U 0 = 0.01 m were chosen the same as in Section 2.1 for the sinusoidal initial functions. In both cases of rectangular and meander-type pulses the number of roots of transcendetal equations (2.24)-(2.26) are the same, N = 8. Dashed vertical lines in each figure show positions of several first roots of corresponding equations. The amplitudes of solitons related to these eigenvalues as per equation (2.23) are presented in Table

Figure 7 .

 7 Figure 7. Cosine initial pulse (dashed line) and the result of its evolution at the distance x = 6 • 10 4 m (solid line).

Figure 8 .

 8 Figure 8. The model amplitude spectrum of random wave field shown in Figure 1.

Figure 9 .

 9 Figure 9. Disintegration of one of the initial random pulses (line 1) on solitons (line 2).

Figure 10 .

 10 Figure 10. Frame (a) -the histogram of soliton distribution (piece-linear line) and the model distribution function (smooth line) for the illustrative example. Frame (b) -the cumulative distribution function versus soliton amplitude corresponding to the model distribution function shown in frame (a). Smooth lines in both frames were obtained by means of spline interpolation of numerical data.

Figure 11 .

 11 Figure 11. Fourier spectrum of wind waves generated in the laboratory tank at two different distances from the ventilator. Solid line 1 represents the spectrum at the distance 100 cm; dashed line 2 represents the spectrum at the distance 300 cm.

Figure 12 .

 12 Figure 12. Fragments of low-frequency components of surface perturbations generated by wind in two spatial points of measurements in the tank. Frame (a) -surface perturbation at the distance 100 cm; frame (b) -at the distance 300 cm.

Figure 13 .

 13 Figure 13. The histogram of soliton numbers versus soliton amplitudes for the time series shown in Figure 12. Line 1 in each frame represent experimental data, and lines 2 -the best fit of these data by the Poisson distribution functions with the parameters λ = 3.85 for the histogram shown in frame (a) and λ = 4.94 for the histogram shown in frame (b).

Figure 14 .

 14 Figure 14. The integral (cumulative) distribution functions (lines) for the experimental data shown in Figure12. Lines 2 represent the best fit approximation with the Poisson cumulative distribution functions with the same parameters λ as in Figure13.

  demonstrates the histogram of soliton numbers versus amplitudes for the time series shown in Figure 12. The experimental data can be approximated by the Poisson distribution function P (n) = λ n e -λ /n!, where the parameter λ = 3.85 for the histogram shown in frame (a) and λ = 4.94 for the histogram shown in frame (b).

Figure 15 .

 15 Figure 15. The fragments of numerical calculations with the input data taken from Figure12illustrating the soliton gas density in both time series.

Acknowledgements. This work was initiated while one of the authors (Y.S.) was the invited visitor at the Laboratoire IRPHE-IOA, Marseille, France several years ago. Y.S. is grateful to CNRS, France for the invitation and to the Laboratory staff for the hospitality. Y.S. also acknowledges the funding of this study from the State task program in the sphere of scientific activity of the Ministry of Education and Science of the Russian Federation (Project No. 5.1246.2017/4.6) and grant of the President of the Russian Federation for state support of leading scientific schools of the Russian Federation (NSH-2685.2018.5). The authors are thankful to Efim Pelinovsky for useful discussions and valuable advices.


Etienne Borré

Jean-François Stumbé

Stéphane Bellemin-Laponnaz

Matteo Mauro
email: mauro@unistra.fr

Control of the photoswitching properties in supramolecular metallopolymeric gels by tuning the coordination metal

Two novel supramolecular metallo-heteropolymers bearing a photo-isomerizable telechelic bis-terpyridine ligand and either Fe(II) or Co(II) coordination metals were synthesized. Both polymers induced gelation of organic solvents at concentration as low as 0.12 wt.% yielding thixotropic gels. Judicious choice of electronic and photophysical properties of both ditopic ligand and metal ion enabled photomechanical response of the supramolecular organogels upon UV light irradiation through molecular design.

The development of supramolecular functional materials represents a challenging research topic of great current interest from both fundamental and application points of view. [1][2] Precise management of weak intermolecular interactions (i.e., supramolecular approaches) in combination with the use of functional building blocks allows not only preparation of self-assembled and hierarchically organized complex systems, but also enables engineering of so-called "smart" materials, such as stimuli-responsive and adaptive soft materials. [3][4][5][6][7] In this respect, we have recently reported on a novel class of supramolecular metallo-polymers based on bisterpyridine ditopic ligands functionalized with photoswitchable diazo-aryl groups. These polymers are capable of inducing solvent gelation affording to photo-responsive and adaptive organo- [8] and hydrogels [9] . The use of zinc(II) salts with terpyridine (tpy) ligands, i.e. a coordination node with relatively high binding constant and kinetic lability, enabled high degree of polymerisation and dynamic metal-ligand bonds thus contributing to self-healing capability of the resulting gels.

Judicious choice of the nature and oxidation state of the coordinating metal centre embedded into the supramolecular architecture is expected to profoundly influence the properties of the resulting metallo-polymers and responsive materials. Indeed, on one hand, one might modify both kinetic lability and thermodynamic stability of the metal(tpy)2 unit (tpy = 2,2':6',2"-terpyridine type of ligands), which are known to control reversibility of the supramolecular self-assembly process between the tpy moiety and the metal cation (i.e., possibility to self-correct and self-repair) as well as solution equilibria, speciation (i.e., mono-tpy vs. bis-tpy species) and degree of supramolecular polymerisation, respectively. [10][11][12][13][14][15][16][17] On the other hand, one may introduce added-value features including photophysical, magnetic and reversible redox properties upon employment of different and properly selected coordination metals. [18][19][20] In particular, 3d transition metals Fe(II) and Co(II) are already known to form stable homoleptic complexes with terpyridine-type of ligands, being logb = 20.9 and 18.3 for Fe(tpy) 2+ and Co(tpy) 2+ , respectively. [21][22] This ensures high degree of polymerisation in the metal-ligand supramolecular self-assembly process under practically-useful concentration conditions. These two transition metals are characterized by a lower kinetic lability than Zn(II) counterparts in the M-tpy interaction, the half-life of M(tpy) 2+ species being in the order Fe(tpy) 2+ > Co(tpy) 2+ > Zn(tpy) 2+ (8400, 50 and 0.1 min in H2O, respectively). [23][24] Aiming at the preparation of supramolecular materials in which photo-mechanical response could be finely controlled by molecular design of constituent building blocks, we investigated the effect of the variation of the metal(II) ion onto the photochemical behaviour of metallopolymeric gels based on photo-active ligand A and B (Scheme 1). Thus, we have chosen Fe(II) and Co(II) on the basis of thermodynamic and kinetic considerations abovementioned as well as on their different ground and excited state electronic properties that are expected to deeply influence the photo-switching ability of the azobenzene moiety encoded into the supramolecular architecture of the metallopolymers. [25] Additionally, the choice of using Fe(II) and Co(II) allows us direct comparison with the previously reported systems based on zinc, given their relatively closer binding constants when compared to 4d metals such as Rh(III) and Ru(II). Furthermore, we decided to use ligand A and B owing to high degree of conjugation between the photoactive moieties (i.e., azoaryl and phenylene-ethynylene, respectively) and the metal-ligand centre to keep highly coupled electronic interactions between them. [START_REF] Burnworth | 27 Detailed speciation studies in solution were previously reported and they can[END_REF] Scheme 1. Schematic synthetic pathway employed for the preparation of the Fe(II) and Co(II) telechelic metallopolymers with ditopic ligand A and B. In the case of polymer Fe-P and Co-P the stoichiometry is m = 2 and n = 1.

The chemical structure of ligand A and B and the schematic synthetic pathway for the preparation of the supramolecular metallo-heteropolymers are displayed in Scheme 1.

Firstly, supramolecular polymerisation process by employing either Fe(BF4)2 or Co(BF4)2 as metal sources and a combination of ligand A and B has been followed in-situ by UVvisible spectroscopy to elucidate the possibility to form metallopolymer and to ascertain the position of the corresponding chemical equilibrium. Such a study was made possible by the peculiar features associated with the formation of each of the M(tpy) 2+ species, where M is either Fe(II) or Co(II). It is important to notice that the BF4 -counter-anion (i.e., a non-coordinating one) has been chosen to favour formation of bis-terpyridine M(tpy) 2+ nodes vs. monomeric M(tpy)X2 (where X = anion or coordinating solvent molecule). The latter acts as chain termination species, thus lowering the degree of polymerisation. [27] UV-visible spectra recorded upon addition of different aliquots of Fe 2+ into a CHCl3:CH3CN 9:1 solvent mixture containing equimolar amount (concentration = 1´10 -5 M) of ditopic ligand A and B are displayed in Figure 1. M in CHCl3:CH3CN 9:1 with Fe(BF4)2 followed by UV-vis spectroscopy (top). The spectra corresponding to 0 and 1 eq. of Fe 2+ added are displayed as black and red trace, respectively. Black arrows indicate the variation of the spectra features upon addition of the salt; plot of the absorbance variation recorded at 600 nm as a function of Fe(BF4)2 added (bottom). Dotted grey lines are only intended to guide the eye.

The absorption spectrum before addition of Fe(II) salt (black trace) is characterized by two intense bands centred at 326 and 370 nm, corresponding to the overlapping singletmanifold ligand-centred ( 1 LC) transitions of both ligand A and B, i.e. 1 p-p* transitions of the terpyridine, phenyleneethynylene and azo moieties. The formally symmetryforbidden weak n-p* transition of the azo moiety is located at ca. 475 nm. Upon addition of Fe 2+ aliquots, the rise of a band at labs = 330 nm is noticed, attributable to the 1 LC band of the p-conjugated ligand coordinated onto Fe 2+ centres. Even more importantly, an intense band in the visible spectrum at 598 nm is observed, which is attributed to the 1 MLCT of the Fe(tpy) 2+ type of complex. [28] Such a band is bathochromically shifted when compared to the Fe(tpy) 2+ benchmark complex as consequence of the more extended p-conjugated system in the tridentate ligand A and B. Interestingly, the spectra display three isosbestic points at labs = 313, 362 and 387 nm suggesting that coordinating moieties are arranged in only one fashion.

Plotting intensity of such 1 MLCT transition over increasing Fe 2+ :(ligand A+B) ratio allowed following the metal-directed supramolecular polymerisation process. The plot clearly show a monotonic increase of the intensity associated to the 1 MLCT band upon increasing Fe 2+ concentration until 1:1 ratio is reached indicating that in such condition complete metal coordination is achieved and the chemical equilibrium is shifted toward the formation of bis-terpyridine species (see Figure 1). As expected for a mechanism involving step-wise coordination of two terpyridine onto a M 2+ ion in which K2 >> K1, [START_REF]tpy + M 2+ ⇄ [M(tpy)] 2+ and tpy + [M(tpy)] 2+ ⇄ M(tpy)2] 2+ coordination event, respectively[END_REF][START_REF] Friese | [END_REF] excess of metal cation did not result in an absorbance variation owing to the higher binding constant of Fe 2+ for the formation of [Fe(tpy)2] 2+ -type species. This demonstrate that the resulted metallopolymer is stable and do not undergo depolymerisation in the presence of excess of metal salts Moreover, the rigid nature of the telechelic ligands A and B points towards formation of linear supramolecular coordination polymers rather than self-assembled metallomacrocycles. [14,31] Titration experiments carried out under identical condition in terms of solvent and ligands concentration, but employing Co(BF4)2 as metal source, afforded similar results, as shown in Figure S1 of the Electronic Supplementary Information (ESI). Upon addition of Co 2+ aliquots, the resulting absorption spectra display appearance of a band at 331 nm attributable to the 1 LC transitions. Noteworthy, three much weaker bands appear at labs ca. 400, 460 and 525 nm that can be ascribed to the MLCT and MC transitions of the Co(II)-ligand complex. [32] Similarly to what observed for the formation of metallosupramolecular coordination polymer with Fe 2+ , and as one can expect on the basis of thermodynamic metal-ligand coordination constants for [Co(tpy)2] 2+ , supramolecular polymerisation was completed at metal:ligand(A+B) ratio of 1 and excess of metal cation did not resulted into formation of dinuclear [Co(ditopic ligand)Co] 4+ species (see Figure S1 of the ESI).

Metallopolymers were synthetized in bulk by refluxing a CHCl3 solution of ligand A:B in a 2:1 ratio and stoichiometric amount of M(BF4)2 (M = Fe 2+ and Co 2+ ) overnight, and the corresponding supramolecular polymers, namely Fe-P and Co-P, were obtained as orange and deep-violet solid, respectively (see Scheme 1, with m = 2 and n = 1). At a first attempt, fixed ratio of 2:1 between ditopic ligands A and B was chosen on the basis of our previous experiments. [8][9] The details of the synthetic procedures are reported as ESI and the corresponding 1 H NMR spectrum for Fe-P displayed in Figure S2 of the ESI. We assume that the formation of metallopolymers results into random organization of both ligand A and B along supramolecular polymers chain. Both Fe-P and Co-P metallopolymers dissolved in DMF were found capable to induce gelation upon addition of either EtOH or CH2Cl2 playing the role as nonsolvent, at overall polymer concentration as low as 0.12 wt.%. This enabled investigation of the effect of electronic properties of the coordination metal on the photo-switching ability of the organogel.

Upon mixing either Fe-P or Co-P metallopolymer at concentration of 20 mg mL -1 in DMF and addition of a 20-fold volume of non-solvent, highly coloured and stable gels were obtained within three hours at room temperature, indicating the high tendency of such supramolecular metallopolyelectrolytes to act as gelating agent though establishment of hydrophobic p-p interactions between p-conjugated systems.

Mechanical properties of both Co-P and Fe-P gels prepared at 0.5 wt.% in DMF:EtOH 1:5 v/v where investigated by means of rheological experiments. As shown in Figure S3, a nonlinear dependency of storage and loss moduli, namely G´ and G´´, over strain was observed for both gels at a constant applied frequency of 1 Hz, with yielding point at about 10% of strain. After such point, further increase of strain gave a sol-gel transition resulting in G´´ exceeding G´ values. Although the presence of more stable and less labile Fe(tpy) 2+ and Co(tpy) 2+ units within the supramolecular polymer skeleton of Co-P and Fe-P when compared to Zn(II) counterparts, [8][9] thixotropic properties of the prepared gels were investigated to elucidate the effect of the interchain p-p interactions onto gel stability. The results are displayed in Figure S4 in which low and high strain values were cycled over time. As shown, while at low (0.5%) strain, storage modulus values exceeded the loss modulus, application of much higher (200%) strain resulted in rupture of the gel network, most likely owing to the breakage of interchain p-p interactions, and subsequent formation of the sol state (G´´ > G´). Interestingly, fast recovery of the gel state could be reversibly obtained when strain was set back to 0.5%, confirming the thixotropic properties of Co-P and Fe-P gels.

In addition, the presence of photo-switching unit A within the metallopolymer structure prompted us to investigate the photochemical behaviour of the so-prepared organogels. Noteworthy, upon irradiation of gel samples with a LED lamp as UV light source (lexc = 365 nm), it was found out that Co-P gel showed neat light-induced contraction over an irradiation time of 12 hours, which is accompanied by syneresis (Figure 2). Indeed, electronic properties of the [Co(tpy)2] 2+ node are such, i.e. weak MLCT and MC band into the visible region, that they can not detrimentally affect the photo-switching of ligand A units.

Figure 2.

Photoresponse under UV irradiation with a continuous-wave LED lamp (lexc = 365 nm, irradiance: 21 mW cm -1 ) of the supramolecular gels containing either Co-P (left) or Fe-P (right). Organogels were prepared in DMF:EtOH (1:20 v /v) at concentration of 0.12 wt.% of metallopolymer as the gelating agent.

On the other hand, no photomechanical response for gel made of Fe-P metallopolymer under identical excitation condition was observed even under an irradiation time as long as 24 hours (see Fig. 2). Such absence of gel response upon light-excitation into the azobenzene p-p* band is ascribed to ultrafast energy transfer (EnT) process from the excited donor azobenzene onto the acceptor MLCT band of the [Fe(tpy)2] 2+ moiety. [33][34] This competitive photophysical process is thus expected to efficiently quench the photo-isomerization of the A unit by funnelling excitation energy into the iron(II)-tpy MLCT state, which in turn deactivates non-radiatively, most likely via 3,5 MC state to the electronic ground state. At the present stage, it is important to notice that the photoswitching ability of the Co(II) metallopolymers enabled surface 2D nano-patterning and monitoring of the photochemical process at the solid-liquid interface by means of scanning probe techniques. [35] In conclusion, two novel metallosupramolecular polymers have been synthetized in a straightforward synthetic pathway by using Fe(II) or Co(II) as coordination metal nodes for the telechelic ditopic ligands. The functional supramolecular metallo-polyelectrolytes, which include photo-isomerizable azoaryl units, are able to induce gelation of organic solvents resulting into thixotropic organogels. Proper tailoring the nature of the coordinating metal enabled the control of the photo-switching ability of the electronically coupled azo moiety. Thus, modulation of photomechanical response in supramolecular gel was achieved in an unprecedented manner. To the best of our knowledge, these represent the first examples of photoswitchable metallopolymeric organogels in which photomechanical response can be controlled by choice of the coordination metal. These results pave the way for better engineering of photo-responsive and adaptive supramolecular materials by chemical design. The University of Strasbourg Institute for Advanced Study (USIAS) is kindly acknowledged for financially supporting this work.

Figure 1 .

 1 Figure 1. Titration of ligand A and B at concentration of 1´10 -5M in CHCl3:CH3CN 9:1 with Fe(BF4)2 followed by UV-vis spectroscopy (top). The spectra corresponding to 0 and 1 eq. of Fe 2+ added are displayed as black and red trace, respectively. Black arrows indicate the variation of the spectra features upon addition of the salt; plot of the absorbance variation recorded at 600 nm as a function of Fe(BF4)2 added (bottom). Dotted grey lines are only intended to guide the eye.


Giulia De Meijere

Vittoria Colizza

Eugenio Valdano

Claudio Castellano

The effect of delayed awareness and fatigue on the efficacy of self-isolation in epidemic control

The isolation of infectious individuals is a key measure of public health for the control of communicable diseases. However, involving a strong perturbation of daily life, it often causes psychosocial distress, and severe financial and social costs. These may act as mechanisms limiting the adoption of the measure in the first place or the adherence throughout its full duration. In addition, difficulty of recognizing mild symptoms or lack of symptoms may impact awareness of the infection and further limit adoption. Here, we study an epidemic model on a network of contacts accounting for limited adherence and delayed awareness to self-isolation, along with fatigue causing overhasty termination. The model allows us to estimate the role of each ingredient and analyze the tradeoff between adherence and duration of self-isolation. We find that the epidemic threshold is very sensitive to an effective compliance that combines the effects of imperfect adherence, delayed awareness and fatigue. If adherence improves for shorter quarantine periods, there exists an optimal duration of isolation, shorter than the infectious period. However, heterogeneities in the connectivity pattern, coupled to a reduced compliance for highly active individuals, may almost completely offset the effectiveness of self-isolation measures on the control of the epidemic.

I. INTRODUCTION

A pillar of non-pharmaceutical interventions for the control of COVID-19 pandemic is the isolation of individuals testing positive for SARS-CoV-2 infection. The aim is to avoid onward propagation of the disease, while contacts are traced to further break the chains of transmission [1]. This measure, however, is met with a set of challenges, as it has no immediate benefit for the index case, but a number of downsides. It often causes psychosocial distress [2], and it may have severe financial and social costs impacting daily life, if a structured support program is not in place.

Ideally, the measure should cover the entire duration of the infectivity period. In practice, isolation may start when a person is already infectious, typically at the onset of symptoms, or when alerted by a contact tracing investigation. Also, the length of the infectious period may be strongly variable across individuals [3,[START_REF]Report of the who-china joint mission on coronavirus disease 2019 (covid-19)[END_REF]. Additional factors may undermine the effectiveness of isolation. Mild symptoms or lack of symptoms may ruin the motivation to respect it, as physical conditions are not an impediment to carry out the daily routine. The legal enforcement of the measure may create tradeoffs discouraging individuals to self-declare as cases [START_REF] Lucas | [END_REF]. Survey data report that adherence is low [6][7][8]. Among the reported reasons for non-adherence are lower socioeconomic grade, psychological distress, inadequate information and long quarantine duration [2].

During COVID-19 pandemic, the duration of isolation has been one flexible component that authorities adapted from initial estimates of 14 days [9] to shorter periods to make the measure more bearable, at the first signs in summer 2020 showing the difficulty of implementation of the measure [START_REF] Silva | U.s. surgeon general blames 'pandemic fatigue' for recent covid-19 surge[END_REF][START_REF] Rahmandad | [END_REF]. Variable durations mark a tradeoff between a long enough period of isolation to prevent onward transmission, and a short enough period that is acceptable by the population. Some countries went as low as 5 to 7 days to increase adherence [START_REF]The Irish Times[END_REF][START_REF] Figaro | [END_REF], especially in countries where self-isolation was not legally compulsory. Further changes (extension to 10 days [START_REF] Lci | [END_REF]) occurred later because of the circulation of the Alpha variant (B.1.1.7 lineage), showing the complexity of the biological and social aspects of setting this public health measure [START_REF] Ashcroft | [END_REF].

As gaps in any of these aspects may undermine the effectiveness of isolation in aiding epidemic control, here we study through mathematical modeling the role of delayed awareness in entering into isolation and fatigue inducing early release of the measure. Our model is a variation of the standard susceptible-infected-susceptible (SIS) compartmental model for infectious disease dynamics [16,[START_REF] Anderson | Infectious Diseases of Humans: Dynamics and Control[END_REF], allowing for three additional compartments: an isolated (Q) compartment, an undecided (U) compartment, and a fatigued (F) compartment. Here we borrow the classical notation Q commonly used in compartmental models to define the isolation of infectious individuals, notably with the susceptible-infected-quarantined-susceptible (SIQS) model [START_REF] Hethcote | [END_REF][19][20][21][22][23]. We do not consider the quarantine as the preventive isolation of suspect cases or of contacts of confirmed cases [24][25][26] , and in the following we will use the terms quarantine and self-isolation as synonyms. The existence of (temporary) immunity against SARS-CoV-2 would suggest the consideration of a SIR-like dynamics. We prefer to consider a SIS-based modeling framework as the absence of an immune state allows us to keep analytical derivations simpler while still providing (a worse case scenario) intuition on the behavioral mechanisms related to the self-isolation measure. We expect that similar results would be obtained for a SIR-like model.

The undecided compartment U corresponds to an intermediate state, following infection, during which awareness arises around the knowledge of being infected, involving a delay before the decision to comply with isolation. This state may also be interpreted as the time between infection and testing (thus corresponding to logistical delays in accessing and performing a test, and obtaining the test results), or to the time between infection and symptoms onset (thus corresponding to a pre-symptomatic phase) [1,27]. Another addition to the standard SIQS model is the possibility that the individual exits selfisolation before its full duration and while still infectious. The compartmental model and transitions are fully explained in the next section.

We investigate the model on a networked population with a mean-field approach, highlighting how the key parameters describing the epidemic dynamics -i.e. the epidemic threshold and the prevalence of infected individuals in the endemic state -depend on the different durations associated to these states. We then consider increasingly more accurate mean-field types of approach, allowing to analyze in detail how individual heterogeneities influence collective properties of the system. We therefore rely on proven effective analytical and numerical tools in order to quantitatively uncover the role of various kinds of imperfections of self-isolation in the spread of a pathogen, which can be of public health relevance to the control of the currently ongoing pandemic.

II. EPIDEMIC MODEL WITH QUARANTINE, DELAY AND FATIGUE

The model we consider is a modification of the usual SIS dynamics [16], based on the existence of three additional compartments, beyond the standard S (susceptible) and I (infected) states.

The contact of a susceptible individual with an infectious one leads to a transmission of the pathogen at rate β. The newly infected individual enters the state U (undecided) preceding the decision on whether to self-isolate or not. The decision process is assumed to be Poissonian with rate µ U . After the decision, the individual enters the quarantined state Q with probability p Q (quarantine probability). With complementary probability he/she instead enters the infected state I. In the latter case the individual behaves as in the standard infected state of the SIS model. In the Q state instead the individual has no contacts, and does not transmit the infection. Compliance to isolation may however end before full recovery, as fatigue sets in. We model this assuming that each quarantined individual transitions to compartment F (fatigued) at rate µ Q . Individuals in state F are infectious with the same transmissibility β of those in state I.

As the progression of the disease does not depend on the isolation status, spontaneous recovery transitions occur from states U, I, Q and F to the susceptible state S, at same rate µ. As a consequence of this choice, the average time spent before reaching the S compartment from any of the infected states U, I, Q and F is 1/µ, also when there exist multiple paths to recovery. The computation of the average time should indeed receive the contributions of all the possible paths, each weighted by the probability of being chosen; when multiple transitions out of a compartment are possible, the average time of each transition must be conditioned on the fact that the other possible transitions were not undergone.

Our model overall depends on five independent parameters p Q , β, µ U , µ Q and µ. The compartments and the transitions allowed between them are depicted in Fig. 1. All transitions are spontaneous except the one taking individuals in state S to the undecided state U, which occurs because of a contact between a susceptible individual and an infectious one. Note that in this model individuals in compartments I, F, U are all infectious with the same transmissibility.

III. MEAN-FIELD APPROACH

Let us define as S(t), U (t), I(t), Q(t), F (t), the probabilities for an individual to be in the respective compartments. The sum of these probabilities equals 1, leaving only 4 independent quantities. We assume a homogeneous pattern of interactions, with average number of contacts ⟨k⟩. Then, the differential equations describing the evolution of the aforementioned probabilities read as follows:

           U = β ⟨k⟩ (1 -I -U -F -Q)(I + U + F ) -(µ U + µ)U İ = µ U (1 -p Q )U -µI Q = µ U p Q U -(µ Q + µ)Q Ḟ = µ Q Q -µF (1)
The disease-free state (S, U, I, Q, F ) = (1, 0, 0, 0, 0) is always an equilibrium solution of the system. Linearization around it shows that it is not stable if λ = β/µ is above the critical value λ c , marking the existence of an endemic state:

λ c = 1 ⟨k⟩ 1 1 -p Q µµ U (µ Q +µ)(µ U +µ) = 1 ⟨k⟩ 1 1 - p Q 1+ T T Q 1+ T U T , (2)
where the temporal scales T = 1/µ, T Q = 1/µ Q and T U = 1/µ U are the average times spent in the corresponding states. Eq. ( 2) contains several known results for limit values of its parameters. p Q = 0 -individuals never in isolation -yields the well-known SIS Mean-field result λ c = 1/ ⟨k⟩. The same limit is recovered if the quarantine has vanishing duration (T /T Q → ∞) or when the time to take a decision diverges (T U /T → ∞).

For generic values of the parameters Eq. ( 2) can be written as

λ c = 1 ⟨k⟩ 1 1 -p ef f Q (3)
where

p ef f Q = p Q 1 + T T Q 1 + T U T . (4)
The quantity p ef f Q is a scaling law turning the effect of compartments U, F into an effective probability to selfisolate in a SIQS model. It is smaller than p Q and reflects the reduction in the efficacy of the quarantine due to undecidedness and fatigue. Note that, even for full compliance with the quarantine prescription (p Q = 1), any value T Q < ∞ or T U > 0 is sufficient to make the threshold finite. For relatively large decision time (compared with recovery time) or small quarantine duration, the factor multiplying p Q in Eq. ( 4) is small and therefore the increase of the epidemic threshold for a full quarantine probability (p Q = 1) compared to none (p Q = 0), might be very limited.

In the endemic state, the densities of individuals in the various compartments are given by

         U * = 1 1+T /T U λ-λc λ I * = (1 -p Q ) 1 1+T U /T λ-λc λ Q * = p Q (1+T U /T )(1+T /T Q ) λ-λc λ F * = p Q (1+T U /T )(1+T Q /T ) λ-λc λ (5)
where the quantity λ-λc λ is the total density of infected individuals 1 -S * = U * + I * + Q * + F * . The total density of infectious individuals instead is

I * tot = I * + U * + F * = 1 ⟨k⟩ 1 λ c - 1 λ . (6)
For any λ, I * tot depends on p Q , T Q and T U only via the value of the epidemic threshold. Equation ( 6) then indicates that, for a given λ, changing parameters in order to increase the epidemic threshold simultaneously reduces the overall prevalence of infectious individuals. Therefore maximizing the epidemic threshold, minimizing I * tot and maximizing S * = λ c /λ are equivalent procedures.

In general, the two parameters p Q and T Q are likely to be dependent, as the perspective of long isolation may discourage people from isolating. We assume that the quarantine probability is a function p Q (T Q ) of its duration. In particular we expect that p Q decreases as T Q increases. This leads to the existence of an optimal quarantine duration T * Q . For small values of T Q adherence to quarantine is high, but its duration is too short, so that people exit from it when they are still infectious. For large values of T Q isolated individuals recover when in isolation, but compliance is low. An optimal tradeoff exists between these two limits. We want to find the optimal duration of self-isolation T Q that minimizes pathogen circulation, i.e. it either minimizes I * tot or maximizes the epidemic threshold. This occurs when the denominator of the epidemic threshold (Eq. ( 2)) attains its minimum, i.e. for

T Q = T * Q such that dp Q dT Q T Q =T * Q = - p Q (T * Q ) T * Q 1 1 + T * Q T . (7)
We make the relationship between p Q and T Q explicit, making the following minimal assumption:

p Q (T Q ) = 1 1 + cT Q /T , (8)
where the parameter c determines how quickly the probability to enter quarantine decays with its duration. This way compliance tends to be perfect (p Q → 1) for extremely short quarantine, while virtually nobody decides to isolate (p Q → 0) if the duration of self-isolation is much longer than the average recovery time.

The optimal quarantine duration, solution of Eq. ( 7), is then

T * Q T = 1 √ c , (9)
corresponding to an optimal probability to quarantine

p * Q = 1/(1 + √ c
) and to the maximum threshold

λ * c = 1 ⟨k⟩ 1 1 - 1 (1+ √ c) 2 1+ T U T . ( 10
)
We find that for the optimal duration of quarantine T * Q the threshold is a decreasing function of the parameter c (see Fig. 2).

IV. HETEROGENEOUS MEAN FIELD APPROACH

We now allow for the more realistic assumptions of individuals to have heterogeneous contact rates. We start investigating this case by means of the Heterogeneous Mean Field (HMF) approximation [16,28], which assumes that the probability of being in a given compartment only depends on the degree k of an individual. Hence the state of the system is described by the set of variables S k (t),

U k (t), I k (t), Q k (t), F k (t),
where k spans the degree values in the network and normalization to 1 holds for any k. This approach is equivalent to assuming that the underlying contact network among individuals is annealed [29], i.e., connections are fully rewired at each time step while preserving the degree of each node. For the standard SIS dynamics on power-law degree-distributed networks with P (k) ∼ k -γ , the HMF approximation provides very accurate results if γ < 5/2 [30] while for larger values of γ it fails for large systems [30,31]. For simplicity we further assume that the network is uncorrelated so that P (k ′ |k) = k ′ P (k ′ )/ ⟨k⟩.

A. Epidemic parameters independent of k

We first consider the case where all individuals behave in the exact same way so that parameters take fixed values.

The HMF equations are easily written down

       Uk = β(1 -I k -U k -Q k -F k )kΘ -(µ U + µ)U k İk = µ U (1 -p Q )U k -µI k Qk = µ U p Q U k -(µ Q + µ)Q k Ḟk = µ Q Q k -µF k ( 11
)
where Θ is the probability that a neighbor of a given node is infectious in an uncorrelated network

Θ = k ′ k ′ P (k ′ ) ⟨k⟩ (I k ′ + U k ′ + F k ′ ). ( 12
)
At stationarity we have

     0 = β(1 -I k -U k -Q k -F k )kΘ -(µ U + µ)U k 0 = µ U (1 -p Q )U k -µI k 0 = µ U p Q U k -(µ Q + µ)Q k 0 = µ Q Q k -µF k (13) whose solution reads          U * k = λkΘ (1+µ U /µ)(1+λkΘ) I * k = µ U µ (1 -p Q )U * k Q * k = µ U p Q µ Q +µ U * k F * k = µ Q µ µ U p Q µ Q +µ U * k ( 14
)
Inserting the stationary values into Eq. ( 12) we find

Θ = k ′ k ′ P (k ′ ) ⟨k⟩ 1 - p Q (1+µ Q /µ)(1+µ/µ U ) (1 + µ U /µ)U * k ′
(15) A nontrivial solution Θ > 0 only exists if the derivative with respect to Θ of the r.h.s. of Eq. ( 15) [where we substitute U * k ′ by its explicit dependence on Θ using Eq. ( 14)] evaluated for Θ = 0, is larger than 1. This condition allows us to determine the epidemic threshold

λ c = 1 1 - p Q 1+ T T Q 1+ T U T ⟨k⟩ ⟨k 2 ⟩ . (16)
We observe that this threshold is simply the HMF threshold for SIS [28] modulated by a factor that takes into account quarantine probability, duration as well as delay.

The effect of topology factorizes. For a homogeneous network Eq. ( 2) is recovered, since k 2 = ⟨k⟩ 2 . We perform numerical checks of these predictions, by simulating SIS dynamics (using a Gillespie optimized algorithm [32]) on networks built according to the uncorrelated configuration model [33]. In this model, we consider an upper cutoff on the degrees -k ∈ [k min = 3, k max = √ N ] -in order to have an uncorrelated network without multiple and self connections. The epidemic threshold is estimated by finding the value of λ = β/µ at which the susceptibility of the system reaches a maximum value [30]. Such susceptibility is computed for the number of infected individuals in the quasistationary regime (the order parameter of the epidemic phase transition). Of course only surviving runs of the dynamics need to be considered. In order to work with the equivalent of surviving runs, we implemented the so-called Quasistationary State method (QS) [30], for which the dynamics never allows the system to enter the healthy absorbing state.

In the following, we consider networks with an exponent γ = 2.25 of the degree distribution and, unless otherwise specified, with a network size N = 10 5 .

In Fig. 3 we plot the epidemic threshold as a function of the system size N for several values of the probability p Q to enter quarantine. We first note an excellent agreement between the theory (dashed lines) and the simulations (symbols), further increasing as N grows. The value of the threshold decreases as a function of size, as a consequence of the diverging second moment at the denominator of Eq. ( 16). A higher quarantine probability leads to an increase of the threshold but for the present choice of parameter values (T = 1, T U = 0.25, T Q = 1), the effect is not dramatic: even a complete participation to quarantine (p Q = 1) implies only a (slightly more than) two-fold increase in the value of the threshold with respect to the p Q = 0 case. In Fig. 4 we show the dependence of the threshold on the duration of quarantine for various values of the probability p Q . We observe that a longer duration of quarantine leads to a larger epidemic threshold, but the effect is sizeable only provided p Q is quite large. In Fig. 5 we show the dependence of the threshold on the duration of quarantine for various delays T U . Here too the threshold increases smoothly with T Q . If the decision time is much smaller than the time to heal the effect becomes relevant also for reasonable values of T Q .

Finally in Fig. 6 we report the dependence of the total density of individuals in the various compartments as a function of λ, showing a fair agreement between theory and numerics. As long as the parameters of the model are finite and strictly positive, each of these densities carries information on the global state of the epidemic, each of them simply being a fraction of the order parameter.

B. Degree-dependent epidemic parameters

Empirical evidence [34] suggests that people having more contacts or being more active tend to be more reluctant in reducing their interactions to prevent contagion. This may be due to the fact that each interrupted contact carries with it a substantial economical, social and/or psychological cost. Real data also suggest that individuals with high activity are more attractive, which may make it more difficult for them to self-isolate given the numerous solicitations they receive from others [23]. It is then quite natural to believe that also adherence to the prescription to self-isolate may be different (and in particular be suppressed) for people with a large number of contacts. In our framework it is possible to model such a realistic element by assuming that the probability to enter quarantine and/or its duration depend on the degree of the node, a proxy of individual activity. In particular it is reasonable to expect both p Q and T Q to decrease with k.

By repeating the calculations already performed in the case with degree-independent parameters we easily find that the epidemic threshold reads

λ c = ⟨k⟩ 1 - p Q (k) 1+ T T Q (k) 1+ T U T k 2 , (17)
where ⟨X(k)⟩ = k P (k)X(k). Hence, depending on how p Q and T Q behave for large k, quarantine may reduce or not the vulnerability of scale-free networks to epidemics. We test this prediction again by performing simulations on networks built according to the uncorrelated configuration model. For reference, we compare with results obtained with degree-independent parameters tuned to have exactly the same average value of the degree-dependent case.

We first check what happens assuming p Q (k) = k min /k, so that compliance is perfect for nodes having minimal connectivity while it becomes very small for large k. In Fig. 7 we report the behavior of the epidemic threshold as a function of T Q for the degree-dependent case and for a degree independent case such that

p Q = ⟨p Q (k)⟩.
It turns out that the threshold is smaller in the degreedependent case and in particular that it grows much more slowly with T Q . The effect of a long self-isolation of less connected individuals is almost completely offset by the little compliance of nodes of large degree.

We then check what happens instead when

T Q (k) = (T -T U )(k max -k)/(k max -k min ).
Complying individuals with few contacts self-isolate until full recovery, whereas individuals with large connectivities spend a vanishing time in quarantine. This degree-dependent scenario is compared in Fig. 8 to the degree-independent one in such a way that in the latter T Q = ⟨T Q (k)⟩. We find that the possibility for a few hubs to undergo shorter isolation . periods than the average individual of the population lowers the epidemic threshold the more the larger the quarantine probability. The linear interpolation between the extreme behaviours of hubs and poorly-connected individuals seems to be slow enough not to completely offset the advantages put forward by the self-isolation prescription.

V. QUENCHED MEAN FIELD APPROACH A. Epidemic parameters independent of k A more refined approach to the model dynamics on networks is provided by the Quenched Mean Field (QMF) approximation [START_REF] Wang | 22nd International Symposium on Reliable Distributed Systems (SRDS'03)[END_REF][START_REF] Van Mieghem | [END_REF][37] (also known as Individual Based mean-field approach), which takes into account the detailed structure of the network encoded in the adjacency matrix A ij . Defining as I i (t), U i (t), Q i (t), F i (t) the probabilities that node i is in state I, U, Q and F respectively, the evolution of the model is described by the set of 4N equations

           Ui = λ(1 -I i -U i -Q i -F i )• • j A ij (I j + U j + F j ) -(µ + µ U )U i İi = µ U (1 -p Q )U i -µI i Qi = µ U p Q U i -(µ + µ Q )Q i Ḟi = µ Q Q i -µF i (18)
By linearizing around the healthy state and imposing the largest eigenvalue of the Jacobian matrix to be equal to 0 one obtains the epidemic threshold

λ c = 1 1 - p Q 1+ T T Q 1+ T U T 1 ρ(A) , (19)
where ρ(A) is the spectral radius of the adjacency matrix (i.e. its largest eigenvalue). This expression is perfectly analogous to the QMF result for standard SIS dynamics (corresponding to Eq. ( 19) for p Q = 0), for which it is well known [16,38] that the QMF threshold is a lower bound of the true threshold. We expect this to be true also for the present modification of the SIS model. Indeed Figs. [START_REF]Report of the who-china joint mission on coronavirus disease 2019 (covid-19)[END_REF] and 5 show that Eq. ( 19) is a tight lower bound, as it is the case for γ = 2.25 < 5/2. For larger values of γ instead additional nontrivial effects [39] make the QMF estimate inaccurate for large networks. We note also that the HMF predictions are slightly closer than QMF to the numerical results. This better performance of HMF with respect to QMF (which occurs also for SIS [40]) is accidental: the additional approximation introduced by HMF partly cancels the error due to the QMF approximation. It is only for larger values of γ > 3 and larger system sizes that QMF reveals its more accurate qualitative behavior.

B. Degree-dependent epidemic parameters

If the parameter p Q depends on k we obtain that the threshold is

λ c = 1 ρ( Ã) . (20)
The quantity à is a modification of the adjacency matrix of the system

à = A   1 -P Q 1 1 + T T Q 1 + T U T   , (21)
where

P Q = diag(p Q (k(i))).
If instead the parameter T Q depends on k we obtain that the threshold is

λ c = 1 ρ( Ã) . (22)
à is now

à = A   1 -p Q 1 1 + T T Q 1 + T U T   , (23)
where

T Q = diag(T Q (k(i)))).

VI. CONCLUSIONS

In this paper we have defined and used an epidemic model to study the role of delayed case detection/infection awareness, compliance to self-isolation, and fatigueinduced early drop-out on the effectiveness of self-isolation as a non-pharmaceutical intervention. We found that adherence to the prescription to self-isolate once infected scaled up the epidemic threshold compared to the simple SIS result, and delay into entering isolation or early release from it resulted in a reduction of effective adherence to self-isolation. If the propensity to enter self-isolation and the time spent isolated decrease with the individual number of contacts (degree), the low adherence of few well socially connected individuals may undermine the effectiveness of the entire non-pharmaceutical measure against the epidemic. This is a key result as it suggests that these phenomena, empirically found [6,41], may strongly limit the impact of isolation programs on the pandemic, unless specific measures are implemented to overcome these barriers [6].

The applicability of this model to real case scenarios would take advantage from being informed by real data for both the parameters describing the COVID-19 dynamics and behavioral parameters. Whereas the former can be obtained by fitting case incidence time-series, the latter rely on a complex interaction between top-down regulations and behavioral adaptations and are hence harder to be inferred from data. The behavioral parameters are therefore explored rather than fit to data. We plan to include data from surveys in a future development of this study. Also additional ingredients must be considered in the model to increase the applicability of the model to more realistic scenarios. First, a more detailed compartmental structure accounting for the different phases of COVID-19 disease progression, to better account for the interplay of different time periods and include asymptomatic and paucisymptomatic states. These may also result in different behaviors, reducing adherence and increasing early drop-outs, compared to symptomatic cases. Our simplified approach has, however, the advantage of being analytically tractable, therefore providing an immediate solution under certain approximation and offering an intuition into the behavior of the system. In this perspective, the SIS model was preferred to a model with immunity. Second, transitions were modeled with Poissonian probability distributions, whereas many of these processes are generally described by broader distributions [42,43] . In this context, we expect that this approximation may impact the recovery process from the state F differently than from the state I.

Further directions can be considered to expand this approach in future work. Here, we assumed that isolation prevents all contacts. In reality, isolation is never 100% effective, due both to behavioral aspects, and hard living constraints (e.g., household crowding [44]). Different degrees of imperfect isolation can be considered in terms of approximations altering the contact pattern. We did not consider in this study the role of quarantine as preventive isolation of suspect cases or contacts of confirmed cases. There is now a large body of literature on the role of contact tracing in combination to isolation and testing in COVID-19 control [START_REF] Ashcroft | [END_REF], and the importance of speeding up this process through digital tools [24]. Beside contact tracing, the introduction of a compartment describing individuals uncertain about their infection status, but still with recommended self-isolation, constitutes an additional component in limiting adherence, as motivation to selfisolate is reduced in absence of symptoms or of a test result confirmation. These processes are likely to be governed by different parameters of quarantine probability and duration. Finally, adherence to self-isolation may be the result of an individual component, explored here, along with a population component defined by a level of awareness and of risk perception that may evolve over time [45], depending on the evolving epidemic context. This may be an important component contributing to the observed relaxation effects after COVID-19 pandemic wave, possibly resulting in case resurgences [46][47][48].

Figure 1 .

 1 Figure 1. Schematic representation of compartments and transition rates of the model .

Figure 2 .

 2 Figure 2. Dependence of the epidemic threshold on the quarantine duration TQ for various values of the parameter c in Eq. (8). Vertical lines indicate the optimal values from Eq. (9). T = 1, TU = 0.2 and ⟨k⟩ = 0.6.

Figure 3 .

 3 Figure 3. Epidemic threshold λc as a function of the system size N for γ = 2.25 and several values of pQ. Dashed lines represent the theoretical prediction, Eq (16), while symbols are the numerical estimates, obtained from the peak of the susceptibility [32]. T = 1, TU = 0.25, TQ = 1.

Figure 4 .

 4 Figure 4. Epidemic threshold λc as a function of TQ for T = 1, TU = 0.2 and various pQ. Symbols are the results of numerical simulations, dashed lines are the predictions of HMF theory, dotted lines are the predictions of QMF theory.

Figure 5 .

 5 Figure 5. Epidemic threshold λc as a function of TQ for T = 1, pQ = 0.7 and various TU . Dashed lines are the predictions of HMF theory and dotted lines are the predictions of QMF theory while symbols are the results of numerical simulations.

Figure 6 .

 6 Figure 6. Density of individuals in the various compartments in the stationary state as a function of λ for T = 1, TU = 0.2, TQ = 0.33 and pQ = 0.7.

Figure 7 .

 7 Figure 7. Epidemic threshold as a function of TQ for T = 1, TU = 0.25. In the heterogeneous case (circles and bottom lines) pQ(k) = kmin/k. In the homogeneous case (squares and top lines) pQ is the same for all nodes: pQ = ⟨pQ(k)⟩ = 0.66. Dashed lines are for the HMF predictions whereas dotted lines are for the QMF ones.

Figure 8 .

 8 Figure 8. Epidemic threshold as a function of pQ for T = 1, TU = 0.25. In the heterogeneous case (circles and bottom lines) TQ(k) = (T -TU )(kmax -k)/(kmax -kmin). In the homogeneous case (squares and top lines) TQ is the same for all nodes TQ = ⟨TQ(k)⟩ = 0.71. Dashed lines are for the HMF predictions whereas dotted lines are for the QMF ones..