Antonio Degasperis
email: antonio.degasperis@roma1.infn.it

Matteo Conforti

Fabio Baronio

Stefan Wabnitz

Effects of nonlinear wave coupling: Accelerated solitons

Introduction

The propagation of nonlinear waves is a broad research subject, which is of long-standing interest to a broad community of researchers, ranging from applied mathematics to phenomenology in various physical contexts. Mathematical modeling of nonlinear wave propagation generally requires the introduction of nonlinear partial differential equations, whose investigation is a very difficult subject in itself. An epochal step in our understanding of nonlinear waves occurred with the discovery that some of such nonlinear models turn out to be integrable. When looking back about sixty years in history, we may recognize that the beginning of modern soliton theory, namely the discovery by Gardner, Green, Kruskal and Miura [1] of the method of solution of the Korteweg-de Vries (KdV) equation, was triggered by a nonlinear lattice wave problem (the celebrated Fermi-Pasta-Ulam problem [START_REF] Fermi | Studies of Nonlinear Problems[END_REF]). Indeed, besides all complicacies of the spectral methods and algebraic structures which form the hard body of the theory, the soliton solutions remain, for very good reasons, the best known and most important character of the story. In short, soliton solutions are solitary non-dispersive waves which, for scalar wave equations, generally look as bell-shaped profiles which may be simply characterized by few parameters such as their amplitude, width and velocity. Popular examples in 1+1 dimensions are the soliton solution u(x, t) = A cosh 2 [p(xvt)]

(1)

of the KdV equation

u t + c 1 u x + c 3 (u xxx + 6uu x ) = 0, (2)
as well as the envelope soliton solution u(x, t) = A e i(kx-ωt) cosh[p(xvt)]

(3) of the nonlinear Schrödinger (NLS) equation

u t = ic 0 u + c 1 u x + ic 2 (u xx + 2|u| 2 u) . ( 4
)
In these two examples, the pulse amplitude A, width 1/p and velocity v are constant parameters.

Besides, the mutual relations among these parameters, namely A = 2p 2 , v = c 1 + 4c 3 p 2 for the KdV equation, and A = p exp(iθ) , v = -c 1 + 2c 2 k (together with ω = -c 0 + c 2 (k 2p 2 )) for the NLS equation, are distinctive of each corresponding nonlinear wave model, and provide a basic prediction for the observational data. Indeed these solitary waves are really simple objects, in that they behave as free particles. Moreover, the integrability character of the nonlinear wave model equations only comes into play in the process of collision between different solitons.

A natural observation which can be made, is that the above discussed soliton behavior is so simple because the field variable is just a scalar. Therefore, in the search for solitons with non trivial motion properties, one may look for multi-component fields, with the expectation that the coupling between the different wave components could produce a far richer soliton phenomenology. To this end, generalizations of known integrable models were introduced rather soon in soliton theory, by letting the field variable u(x, t) take values in a vector space or in a matrix algebra. However, in these first attempts the soliton solutions were, once again, relatively simple solitary waves. For instance, the matrix version [START_REF] Degasperis | [END_REF] of the KdV equation (2) (here {A, B} is the anticommutator AB + BA)

U t + c 1 U x + c 3 (U xxx + 3{U, U x }) = 0 , (5)
U being a N × N matrix, has the soliton solution

U (x, t) = A cosh 2 [p(x -vt)] P ( 6
)
where P is a constant projector, P 2 = P , and the parameters A, p and v are related to each other as in the scalar case (1). Similarly, in the two component Manakov [4] generalization of the NLS equation u

(1) t = ic 0 u (1) + c 1 u

(1) ) xx + 2(|u (1) | 2 + |u (2) | 2 )u (1) ] u

x + ic 2 [u (1
(2) t = ic 0 u (2) + c 1 u

(2) x + ic 2 [u (2)
xx + 2(|u (1) | 2 + |u (2) | 2 )u (2) ] [START_REF] Calogero | [END_REF] the soliton behavior is the same as in the scalar case (3), namely u (1) (x, t) u (2) (x, t)

= A e i(kx-ωt) cosh[p(xvt)] α 1 α 2 (8) where (α 1 , α 2 ) is a constant unit complex vector,

|α 1 | 2 + |α 2 | 2 = 1.
The decisive step was taken in 1976 by F. Calogero and one the authors (A.D.) [5], who realized that a full generalization of the scalar case requires that also the constant coefficients c 0 , c 1 , c 2 etc. in the evolution equation, see [START_REF] Fermi | Studies of Nonlinear Problems[END_REF] and (4), should be matrices rather than scalars. The first instances of such generalized models were constructed in the KdV-type hierarchy where even the lowest member, which features only first-order derivatives, is a non-trivial one. This is the case of the matrix equation

U t + i[C 0 , U] + 1 2 {C 1 , U x } + [W, U ] = 0 , W x = 1 2 [C 1 , U] , W (-∞, t) = 0 , ( 9
)
where now the constant coefficients C 0 and C 1 (as well as the dependent variables U and W ) are N × N matrices. Moreover the requirement that these coefficients are matrices forces the introduction of the "auxiliary" field W (x, t). The soliton solution of equation ( 9) takes the form

U (x, t) = A cosh 2 {p[x -ξ(t)]} P (t) (10)
where the amplitude A and the width 1/p are the same as before, A = 2p 2 , but the soliton position ξ(t) is no longer a linear function of time as in the case of a free particle. This happens because the soliton speed dξ(t)/dt is related to the projection matrix P (t), which projects on a one-dimensional subspace (hence it is also termed "polarization" matrix) via the equation

dξ(t) dt = Tr[C 1 P (t)] , (11)
while the polarization matrix P (t) itself satisfies the nonlinear ordinary differential equation [START_REF] Degasperis | Nonlinear evolution equations solvable by the spectral transform[END_REF]

dP (t) dt = [pC 1 -iC 0 , P (t)] + 2pP (t)C 1 [1 -P (t)] . (12)
Equation ( 10) describes an accelerated soliton, as the acceleration d 2 ξ(t)/dt 2 is not vanishing, and its motion is that of a particle in a non constant potential. Indeed the simplest instance of equation of the type ( 9) is provided by the so-called "Boomeron equation" [START_REF] Calogero | [END_REF], which in the 2 × 2 case yields the following system describing the coupling of a scalar field z(x, t) with two 3-dimensional vector fields u(x, t) and w(x, t)

z t = b • u x , u t = a ∧ u + z x b + w ∧ u , w x = -2b ∧ u (13)
whose soliton solution reads

z = A/ cosh 2 {p[x -ξ(t)]} , u = z n(t) , w = -2p(1 + tanh 2 {p[x - ξ(t)]}) b ∧ n(t).
Here a and b are arbitrary constant 3-dimensional vectors, and the unit 3dimensional vector n(t) is the polarization of the soliton whose time-dependence is given by the equation

dn dt = a ∧ n -2p(b ∧ n) ∧ n , n • n = 1 , ( 14
)
which, in the special case a ∧ b = 0, coincides with the Landau-Lifshitz equation for ferromagnetic materials. The explicit solution of the evolution equation ( 14) for the polarization vector n(t), together with the expression

dξ(t) dt = -b • n(t) , ( 15
)
of the soliton velocity, shows that the soliton moves as a particle in a potential barrier in such a way that its velocity at very large times becomes constant (as for a free particle) but with different velocities in the past (t = -∞) and in the future (t = +∞). This implies that, in appropriate systems of reference, the soliton in the future goes back to the same region where it came from in the past, like a boomerang. This peculiar behavior motivates the introduction of the term "boomeron" for describing such a soliton. Moreover, and less generically, by appropriately choosing the soliton parameters and the coefficient vectors a and b, the soliton motion may even coincide with that of a particle trapped in a binding potential, with the end result that the soliton position oscillates with periodic motion in a confined space interval. In this particular case, the soliton has been termed "trappon". These findings were initially considered as a mere curiosity, because of the lack of any apparent connection to a physically observable model or application. Quite recently [8], a reduction of the boomeron equation ( 13), the so called "zoomeron" (or Calapso) equation, was found to be related to isothermal surfaces, and indeed it was already introduced almost a century ago in differential geometry by the italian geometer Pasquale Calapso. In the last few years, different integrable wave equations of boomeronic type were constructed and applied to systems of coupled equations, in particular of NLS-type [9,10]. Moreover, the search of models of physical and applicative interest has progressed until several of such models have been found. One such example is the following model (which may be considered as an extension of the Manakov system, see( 7)):

u (1) t = ic 0 u (2) + c 1 u (1)
x + wu (2)

+ ic 2 [u (1) xx + 2(|u (1) | 2 + |u (2) | 2 )u (1) ] u (2) t = ic 0 u (1) -c 1 u (2)
xw * u (1)

+ ic 2 [u (2) xx + 2(|u (1) | 2 + |u (2) | 2 )u (2) ] w x = 2c 1 u (1) u (2) * , w(-∞, t) = 0 , ( 16
)
whose soliton solution (compare with (8))

u (1) (x, t) u (2) (x, t) = A e i(kx-ωt) cosh{p[x -ξ(t)]} α 1 (t) α 2 (t) , |α 1 (t)| 2 + |α 2 (t)| 2 = 1 , ( 17
)
does indeed feature a boomeronic behavior, as implied by the evolution equation of its 2dimensional complex unit polarization vector, namely α(t) = (α 1 (t) , α 2 (t)),

dα 1 dt = ic 0 α 2 + 2pc 1 |α 2 | 2 α 1 , dα 2 dt = ic 0 α 1 + 2pc 1 |α 1 | 2 α 2 (18)
as well as the evolution equation of the soliton center

dξ(t) dt = c 1 (|α 2 (t)| 2 -|α 1 (t)| 2 ) + 2c 2 k , ( 19
)
where, as before, A = p exp(iθ), ω = c 2 (k 2 -p 2 ). Again here the soliton speed is time-dependent since it is related by (19) to the polarization vector (α 1 (t) , α 2 (t)), which satisfies the nonlinear equations (18). Alternatively, the nonlinear ODEs (18-19) may be either directly integrated (with some effort) or, more easily, the solution (α 1 (t) , α 2 (t)) may be obtained from the Lax pair by the Darboux-dressing technique [11,12]. Once again, the soliton (17) features both boomeronic and trapponic behaviors and, more importantly, the system of nonlinear PDEs ( 16) is per se of great interest in nonlinear optics applications. In fact, equations ( 16) model wave propagation with quadratic and cubic Kerr-type coupling, in conjunction with Schrödingertype dispersion. Quite interestingly, a reduced case of the system (16), which is obtained by letting the dispersion coefficient vanish, i.e., c 2 = 0, and by giving the independent variables x and t the exchanged role of time and, respectively, space, is the celebrated integrable model of the interaction of three resonant waves [13,14]. Thus we may conclude that also this well-known model exhibits soliton solutions of boomeronic type [15]. Because of the high physical relevance of the three-wave model, the main focus of this article will be the discussion of the properties of its boomeron-type soliton solutions with particular reference to nonlinear optics. In the next Section 2 we briefly review the general mathematical setting of the construction of matrix integrable wave equations based on the Lax pair approach. We shall limit ourselves to display only one model of coupled wave equations, namely the system (16), which is derived by means of a reduction from general matrix equations. Next in Section 3 we will detail the properties of a family of boomeronic solutions of the three-wave resonant interaction equations. Namely, we will discuss the stability, the asymptotic states and the collisions of three-wave boomerons.

Finally, an outlook for further research and applications of accelerated solitons is presented in Section 4.

Matrix approach to integrable wave coupling

The universal applicability of the NLS equation (4) as a model for the amplitude modulation of waves with weak dispersion and nonlinearity justifies our search for multi-component versions of such equation. The natural approach for constructing integrable multi-component NLS-type equations involves starting with the Lax pair that is associated with the scalar NLS equation (4), and replacing scalar dynamical variables with (generally rectangular) matrices [9]. In fact, we set up now this method in such a way that the scalar NLS equation ( 4) is recovered as just the lowest dimensional case. All systems of coupled PDEs that we may construct in this way [9,10] are special (reduced) cases of the following "mother" matrix PDE:

Q t = [C (0) , Q] + σ[C (1) , Q x ] -σ{Q, W } -ic 2 σ(Q xx -2Q 3 ) , W x = [C (1) , Q 2 ] , ( 20
)
where the dependent variables Q = Q(x, t) and W = W (x, t) are (N (+) +N (-) )×(N (+) +N (-) ) block matrices of the form

Q = 0 N (+) ×N (+) Q (+) Q (-) 0 N (-) ×N (-) , W = W (+) 0 N (+) ×N (-) 0 N (-) ×N (+) W (-) . ( 21
)
In self-evident notation, the diagonal entries W (+) , W (-) are square matrices of dimension N (+) × N (+) and N (-) × N (-) , while the two off-diagonal rectangular blocks Q (+) and Q (-) , are N (+) × N (-) or N (-) × N (+) matrices, respectively, where N (+) and N (-) are arbitrary, positive integers. Here, and in the following, the constant diagonal matrix σ in ( 20) is, in self-evident notation,

σ = 1 N (+) ×N (+) 0 N (+) ×N (-) 0 N (-) ×N (+) -1 N (-) ×N (-) , (22)
while C (0) , C (1) are arbitrary constant block-diagonal matrices,

C (j) = C (j)(+) 0 N (+) ×N (-) 0 N (-) ×N (+) C (j)(-) , j = 0, 1, (23)
and the block-diagonal matrix W = W (x, t) is an auxiliary dependent variable. As usual [A, B] and {A, B} are the commutator AB -BA and, respectively, the anticommutator AB + BA.

The constant scalar coefficient c 2 is a real dispersion parameter.

The matrix equation (20) represents the compatibility condition for the Lax pair

ψ x = Xψ , ψ t = T ψ , ( 24
)
where ψ, X and T are (N (+) + N (-) ) × (N (+) + N (-) ) square matrices, where ψ = ψ(x, t, k) is a common solution of the two linear ordinary differential matrix equations (24), while X = X(x, t, k) and T = T (x, t, k) depend on the coordinate x, the time t and the complex spectral parameter k according to the definitions

X(x, t, k) = ikσ + Q(x, t) , ( 25
)
T (x, t, k) = 2c 2 k [ikσ + Q(x, t)] + 2ikC (1) + ic 2 σ Q 2 (x, t) -Q x (x, t) -σW + +σ[C (1) , Q] + C (0) . ( 26
)
In order to simplify the notation, from now on (and as already done in the previous formula) we do not specify the dimension of the matrices 0 and 1, and we may even omit to write the matrix 1 altogether, as we trust the reader will not be confused by this omission.

Let us now consider of the condition

Q † (x, t) = S Q(x, t) S ( 27
)
on the solution Q(x, t) of the matrix evolution equation (20) where the dagger stands for hermitian conjugation. The constant matrix S is block diagonal,

S = S (+) 0 0 S (-) ; (28)
its off-diagonal blocks are vanishing rectangular matrices, while its diagonal blocks S (+) and S (-) are N (+) × N (+) or N (-) × N (-) , respectively, diagonal matrices whose diagonal elements s

(±) n
(with no loss of generality) are signs. Namely,

S (±) = diag (s (±) 1 , • • • , s (±) N (±) ) , s (±) n 2 = 1 . ( 29
)
This of course implies the relations S 2 = 1, S (+) 2 = 1, S (-) 2 = 1. The reduction equation ( 27) is well motivated by the fact that it captures several interesting models of dispersive propagation of multicomponent waves in weakly nonlinear media. To this aim, it is convenient to rewrite the matrix PDEs (20) in terms of the blocks Q (+) , Q (-) , W (+) and W (-) , see (21). These read

Q (±) t = C (0)(±) Q (±) -Q (±) C (0)(∓) ± C (1)(±) Q (±) x -Q (±) x C (1)(∓) ∓ W (±) Q (±) + Q (±) W (∓) ∓ ic 2 Q (±) xx -2 Q (±) Q (∓) Q (±) , ( 30
)
W (±) x = C (1)(±) , Q (±) Q (∓) . ( 31
)
where C (j)(+) and C (j) (-) , are the N (+) × N (+) , or N (-) × N (-) , respectively, constant square matrix blocks of C (j) , see (23). The reduction condition ( 27) is accounted for by introducing the dependent variable U (x, t) through the definitions

Q (-) (x, t) = U (x, t) , Q (+) (x, t) = S (+) U † (x, t)S (-) . ( 32
)
Since the constant c 2 is real, the above expressions of Q (+) and Q (-) in terms of the single variable U (x, t) are compatible with the equations (30), which then reduce to

U t = C (0)(-) U -U C (0)(+) -C (1)(-) U x -U x C (1)(+) + + W (-) U + U W (+) + ic 2 U xx -2 U S (+) U † S (-) U , ( 33
)
W (+) x = C (1)(+) , S (+) U † S (-) U , ( 34
)
W (-) x = C (1)(-) , U S (+) U † S (-) . ( 35
)
Here U is an N (-) ×N (+) rectangular matrix, (see (32) and ( 21)), while the variables W (±) (x, t) are square matrices: W (+) is an N (+) × N (+) matrix and W (-) is an N (-) × N (-) matrix, respectively, and it is easily seen that they satisfy the "hermitianity" conditions

W (+) = -S (+) W (+) † S (+) , W (-) = -S (-) W (-) † S (-) . ( 36
)
Similarly, the constant matrices C (j)(±) satisfy the following conditions:

C (j)(±) = -(-1) j S (±) C (j)(±) † S (±) , j = 0, 1 . ( 37
)
As for the "sign" matrices S (+) and S (-) , see (29), we note that one could set, for instance, s (+) 1

= 1 with no loss of generality, but we prefer to keep the symmetrical, though redundant, notation (29). We observe that the general matrix nonlinear evolution equation ( 33) may be well specialized to quite a large family of coupled-wave equations by playing with various choices of the integers N (+) and N (-) , and of the constant matrix coefficients C (j)(+) and C (j) (-) . This exercise is certainly worth doing since it happens that among this family of wave equations there are some models which look interesting in applicative contexts. While we refer the reader to [9] and [10] for the derivation of other model wave equations, we specify here the values of the integers N (+) , N (-) and the matrix coefficients C (j)(+) and C (j)(-) which correspond to the system (16). This system of two coupled NLS equations is obtained for N (+) = 1 and N (-) = 2, together with S (+) = 1 and S (-) =diag(-1, -1). Here the choice of the matrix coefficients is

C (0)(+) = C (1)(+) = 0, C (0)(-) = 0 ic 0 -ic 0 0 , C (1)(-) = -c 1 0 0 c 1 , ( 38
)
while the dependent variables are the two components (u (1) (x, t) , u (2) (x, t)) of the 2-dimensional (column) vector U and the function w(x, t) is defined through the off-diagonal matrix W (-)

W (-) = 0 w -w * 0 ( 39
)
while W (+) = 0. The case in which the dispersion coefficient vanishes, c 2 = 0, coincides with the model of resonant interaction of three waves. The boomeron solutions of this system can be obtained by applying the Darboux-dressing method [11] to the Lax pair (24): the rich phenomenology of these accelerated solitons is the content of the next section.

Propagation, stability and scattering of parametric three-wave solitons

Three-wave resonant interaction (TWRI) is ubiquitous in various branches of science, as it describes the mixing of waves with different frequencies in weakly nonlinear and dispersive media. The TWRI model is typically encountered in the description of any conservative nonlinear medium where the nonlinear dynamics can be considered as a perturbation of the linear wave solution, the lowest-order nonlinearity is quadratic in the field amplitudes and the phasematching (or resonance) condition is satisfied. After the first milestone studies on TWRI [13,[16][17][18][19][20][21][22], many theoretical investigations were made in the context of nonlinear optics, plasma physics, acoustics, fluid dynamics, and solid state physics. TWRI has been extensively studied in the context of nonlinear optics, since it applies to parametric amplification, second harmonic generation and frequency conversion, stimulated Raman and Brillouin scattering and light speed control [23][24][25][26][27]. TWRI also occurs in plasmas in the saturation of parametric decay instabilities, radio frequency heating and laser-plasma interactions [17,18,[28][29][30]. In the TWRI scenarios, parametric TWRI solitons play a pivotal role because of their particle-like behavior which makes them important mechanisms of energy transport, waves manipulation and processing. The particle-like behavior of solitons is responsible for numerous phenomena, which cannot occur with linear waves. Two classes of analytical soliton solutions of the TWRI have been known for over three decades. The first type of solitons describes the mixing of three pulses which travel with their respective linear group velocity, and interact for just a short time [13,[20][21][22]31]. The second type of solitons, also known as simultons, are formed as a result of the mutual trapping of pulse envelopes at the three different frequencies. Hence the three wave packets travel locked together with a common group velocity [16][17][18][19]. In this section we present a new three-parameter family of bright-bright-dark simultons that travel with a common, locked velocity [32]. The most remarkable physical property of the present simultons is that their speed can be continuously varied by means of adjusting the energy of the two bright pulses. We have investigated the propagation stability of TWRI simultons (TWRIS) and we have found that a stable triplet loses its stability as soon as its velocity decreases below a well defined critical value. Another feature of a TWRIS is that an unstable triplet decays into a stable one through the emission of a pulse, followed by acceleration up to reaching a constant velocity. Moreover, a stable triplet may be excited into an unstable simulton by slowing down as a result of the absorption of an isolated wave [33]. Such processes are exactly described in terms of an analytical higher-order soliton solution with varying speed, or boomeron. The present TWRIS scattering process may be pictured as the interaction of radiation with a twolevel atomic system: transitions among excited and ground soliton states are induced by the absorption and spontaneous emission of a wave. To conclude our investigations on TWRI solitons we have also performed numerical studies on TWRIS collisions, finding that TWRIS with different speeds are stable upon collision, and TWRIS with equal speeds interact with each other in a way which is strongly dependent upon their initial relative phase.

Three-wave interaction equations

For our convenience, let us first rewrite the coupled partial differential equations (PDEs) (16) with vanishing dispersion, c 2 = 0, as [13]:

E 1t -V 1 E 1z = E * 2 E * 3 , E 2t -V 2 E 2z = -E * 1 E * 3 , ( 40
) E 3t -V 3 E 3z = E * 1 E * 2 ,
where the subscripts t and z denote derivatives in the longitudinal and transverse dimension, respectively. Moreover, E n = E n (z, t) are the complex amplitudes of the three waves, V n are their linear velocities, and n = 1, 2, 3. For the reader's benefit, we provide here the transformation of the dependent and independent variables which maps the system (16) into the form (40). The coordinates x and t in ( 16) are related to the new coordinates z and t by the formula

z = 1 2 (V 1 V 3 + V 2 V 3 -2V 1 V 2 )t + 1 2c 1 V 3 (V 2 -V 1 )x , t = 1 2c 1 (V 1 -V 2 )x + 1 2 (V 1 + V 2 -2V 3 )t , (41)
while the relation between amplitudes reads

E 1 = V 2 -V 3 ρ u (1) * , E 2 = V 1 -V 3 ρ u (2) , E 3 = V 1 -V 2 4c 2 1 ρ (w + ic 0 ) , ρ = (V 1 -V 2 )(V 1 -V 3 )(V 2 -V 3 )/(4c 2 1 ) . ( 42
)
Moreover, for notational convenience, we replace in (40) the new time variable t with t. We assume here the ordering V 1 > V 2 > V 3 which, together with the above choice of signs before the quadratic terms, entails the non-explosive character of the interaction. In the following, with no loss of generality, we shall write the equations (40) in a reference frame such that V 3 = 0. A remarkable property of the equations ( 40) is their invariance with respect to the transformation

Ên (z, t) = s exp[i(q n z n + α n )]E n (sz + z 0 , st + t 0 ) (43)
where α 1 + α 2 + α 3 = 0, q n = q(V n+1 -V n+2 ), z n = z + V n t are the characteristic coordinates and n = 1, 2, 3 mod(3). As the transformation (43) depends on six real parameters, namely α 1 , α 2 , s, q, z 0 and t 0 , clearly one may introduce these parameters in the expression of any given solution of the TWRI equation. The evolution equations ( 40) represent an infinite-dimensional Hamiltonian dynamical system, with the conserved Hamiltonian

H = 1 4i +∞ -∞ V 1 (E * 1z E 1 -E 1z E * 1 ) -V 2 (E * 2z E 2 -E 2z E * 2 ) + + V 3 (E * 3z E 3 -E 3z E * 3 ) + 2E 1 E 2 E 3 -2E * 1 E * 2 E * 3 dz, ( 44
)
energies (Manley-Rowe invariants)

I 12 = I 1 + I 2 = 1 2 +∞ -∞ (|E 1 | 2 + |E 2 | 2 )dz, ( 45
)
I 23 = I 2 + I 3 = 1 2 +∞ -∞ (|E 2 | 2 + |E 3 | 2 )dz, (46)
and total transverse momentum

J = 1 4i +∞ -∞ [(E * 1 E 1z -E 1 E * 1z ) -(E * 2 E 2z -E 2 E * 2z ) + (E * 3 E 3z -E 3 E * 3z )]dz. ( 47
)
Each of the above conserved quantities ( 45)-( 47) is related to a given internal parameter of the TWRIS which, in turn, is associated with a symmetry (e.g., phase rotation or space translation) of the TWRI equations ( 40) [34]. As a consequence, one may expect that equation ( 40) possess a three-parameter family of soliton solutions.

Three-parameters simulton family

In this section, we will discuss the recently found simulton family of solutions of equations ( 40) [32], by using the general results on TWRI equations which were presented in Ref. [15]. The expression of these simultons is

E 1 = 2p a |b| 2 + a 2 g 1 g(V 1 -V 2 ) exp[i(q 1 z 1 -χz + ωt)] cosh[B(z + V t)] , ( 48
)
E 2 = -2p b |b| 2 + a 2 g 2 g(V 1 -V 2 ) exp[i(q 2 z 2 + χz -ωt)] cosh[B(z + V t)] , (49)
E 3 = {1 + 2p b * |b| 2 + a 2 [1 -tanh[B(z + V t)]]} a g 3 exp(iq 3 z 3 ) g(V 1 -V 2 ) (50)
where

b = (Q -1)(p + ik/Q), r = p 2 -k 2 -a 2 , Q = 1 p 1 2 [ r + r 2 + 4k 2 p 2 ], B = p[ V 1 + V 2 -Q(V 1 -V 2 ) ]/(V 1 -V 2 ), V = 2V 1 V 2 /[ V 1 + V 2 -Q(V 1 -V 2 ) ], χ = k[ V 1 + V 2 -(V 1 -V 2 )/Q ]/(V 1 -V 2 ), ω = -2kV 1 V 2 /(V 1 -V 2 ), q n = q(V n+1 -V n+2 ), g n = |(V n -V n+1 ) (V n -V n+2 )| -1/2 , g = g 1 g 2 g 3 , (51)
and n = 1, 2, 3 mod(3). The above expressions depend on the six real parameters V 1 , V 2 , p, k, q, a > 0. From the definition of Q, one can see that the above parameters must be chosen so that if k = 0, then p 2 > a 2 . The TWRIS is composed of two bright pulses (48,49), and a kink or shock-like pulse (50), which travel with a common locked velocity V which does not coincide with any of the three characteristic velocities V 1 , V 2 , V 3 . The expressions (48)-(50) may be represented in a more convenient form as

E n (ξ, τ) = U n (ξ) exp[ iΦ n (ξ, τ) ], n = 1, 2, 3. (52)
Here we use a reference frame which moves along with the soliton, with coordinates ξ = z + V t, τ = t where U and Φ are real functions and Φ n (ξ, τ) = φ n τ +f n (ξ). A simple analysis of (52) shows that, for any value of the parameters, the pulse amplitudes U 1 (ξ), U 2 (ξ) and U 3 (ξ) are even functions of ξ and the phase constants satisfy φ 1 +φ 2 +φ 3 = 0. On the other hand, if k = 0, the phase profiles are all piecewise linear in ξ, and obey the condition f 1 (ξ) + f 2 (ξ) + f 3 (ξ) = 0 for pξ > 0 and f 1 (ξ) + f 2 (ξ) + f 3 (ξ) = π for pξ < 0. Whereas for k = 0 the phase profile f 3 (ξ) is nonlinear and cos[f 1 (ξ) + f 2 (ξ) + f 3 (ξ)] is an odd function of ξ; moreover the kink pulse E 3 is "grey" if k = 0 and is "dark" if k = 0. Such amplitude and phase front profiles prevent a net energy exchange among the three waves. It is important to point out that the condition -1 < Q < 1 leads to a speed V that lies in-between the characteristic velocities V 1 and V 2 of the two bright pulses, i.e. V 1 > V > V 2 . The above described properties mean that TWRIS represent a significant generalization with respect to previously known three-wave solitons which exhibit a simple (constant) phase profile and correspond to the special case k = q = 0 , r > 0 [16][17][18]. In Fig. 1 we plot two examples of TWRIS amplitude and phase-fronts (48)-(50). It is interesting to consider the physical meaning of the various TWRIS parameters appearing in (48)-(50). For a given choice of the characteristic linear velocities V 1 and V 2 , we are left with the four independent parameters p, k, q, and a. We may note that p is basically associated with the scaling of the wave amplitudes, as well as of the coordinates z and t. The parameter a determines the amplitude of the asymptotic plateau of the kink E 3 . The value of k is related to the wave-number of the soliton. The parameter q adds a phase contribution which is linear in z and t. Since the system (40) is invariant under a transformation (43), without loss of generality we may set a = 1, which reduces the number of essential parameters to just three, corresponding to the three symmetries of equations (40). The parameters p, k, q in (48)-(50) may be more conveniently mapped into the parameters V, φ 1 , φ 2 of equation ( 52), which provide a more direct physical insight into the features of a TWRIS. Such a mapping is obtained by comparing equations (48)-( 50) with (52), and reads as:

V = 2V 1 V 2 V 1 + V 2 -Q(V 1 -V 2 ) ,
φ 1 = qV 2 (V 1 -V ) + χV + ω, φ 2 = qV 1 (V -V 2 ) -χV -ω, ( 53
)
whereas the inverse mapping:

q = φ 1 + φ 2 V (V 1 -V 2 ) , k = φ 1 -qV 2 (V 1 -V ) V V1+V2 V1-V2 -V/Q -2V1V2 V1-V2 , p = sign(Q) Q 2 (k 2 + |a| 2 ) -k 2 Q 2 (1 -Q 2 ) . ( 54
)
In Fig. 2 we plot the mapping of the parameters (V, φ 1 ), in the subspace φ 1 = -φ 2 , into the new parameters (p, k) (q = 0, a = 1). It is worth noting that the condition

k 2 > Q 2 a 2 /(1 -Q 2 )
must be satisfied in order to obtain real values of p. This implies that, for a fixed V , there exists a circular boundary for the allowed domain of phase constants φ 1,2 (i.e., the region inside the shaded areas in Fig. 2). Therefore a TWRIS may be simply expressed in terms of its velocity V and the two phase constants φ 1 and φ 2 . Let us investigate what are the TWRIS properties for a fixed choice of the linear velocities V 1 and V 2 , upon variations of its energy flows and transverse momentum. As an example, Fig. 3(a) shows the dependence of the phase constant φ 1 on the locked velocity V , for the case where φ 2 = -φ 1 , with different values of the conserved energy I 12 . Fig. 3(b) illustrates the dependence of the energies I 1 and I 2 (which happen to be timeindependent for a TWRIS) on the locked velocity V , for different choices of the phase constant φ 1 = -φ 2 . Figs. 3(c)-(d) show the dependence of the energies I 12 and of the renormalized transverse momentum J [as defined by (57), see discussion below] upon the locked velocity V respectively, for different choices of the phase constant φ 1 = -φ 2 . As it can be seen from the above figures, the intensity, the phase profiles, as well as the energy distribution among the different wave packets strongly depend upon the value of the locked velocity V .

Three-wave simulton stability

The next crucial issue is the propagation stability of TWRIS. A first insight into this problem may be provided by performing a linear stability analysis as it was done in Ref. [35]. Let us consider a perturbed TWRIS of the form Ẽn (ξ, τ) = (x n (ξ) + P n (ξ, τ))e iφnτ , n = 1, 2, 3 where

x n (ξ) = U n (ξ) exp[if n (ξ)
] is the soliton profile; for a weak perturbation, we take |P n | << |x n |. By inserting the above ansatz in equations (40), and by retaining only linear terms in P n , one obtains a linear system of PDEs. For the numerical analysis, these PDEs can be reduced to a system of ordinary differential equations dP (τ )/dτ = M P (τ ), by approximating the spatial derivatives with finite differences, where P is the perturbation vector sampled on a finite grid. A sufficient condition for the instability of a stationary solution x n (ξ) is that the matrix M has at least one eigenvalue with positive real part. Numerical computations over a wide parameter range show that eigenvalues of M exist with a positive real part whenever p < 0. On the other hand, for p > 0 the largest real part of the eigenvalues is equal to zero, which means that the TWRIS are only neutrally stable. Note that the instability condition p < 0 leads to the inequality

V < V cr = 2V 1 V 2 /(V 1 + V 2 )
. Extensive numerical integrations of equations (40) confirm that TWRIS with V < V cr (V > V cr ) are always unstable (stable). The propagation of either stable or unstable TWRIS is illustrated in Fig. 4, which shows the general feature of unstable solitons with V < V cr . Namely, the simulton decays into a stable soliton with V > V cr , and it emits a pulse in the wave E 3 . It is quite remarkable that the dynamics of the decay from unstable into stable solitons may be exactly described by analytical solutions with variable velocity or boomerons, as it will be shown in the next section. We performed further investigations of TWRIS stability by carrying out a multi-scale asymptotic analysis (MAA) [36][37][38][39][40]. In this way, one obtains the following condition which defines the borderline between stable and unstable TWRIS

G = ∂(I 12 , I 23 , J) ∂(φ 1 , φ 2 , V ) = 0, ( 55
)
where G is the Jacobian of the constants of motion I 12 , I 23 , J with respect to φ 1 , φ 3 , V (φ 3 = -φ 1φ 2 ). Note that in (55), I 23 and J are obtained by re-normalizing the divergent integrals [START_REF] Hasegawa | Solitons in Optical Communications[END_REF] and (47) according to the prescription [41]

I 23 = 1 2 +∞ -∞ (|E 2 | 2 + |E 3 | 2 -|E 30 | 2 )dz, (56)
J = J - 1 4i +∞ -∞ (E * 3 E 3z -E 3 E * 3z ) |E 30 | 2 |E 3 | 2 dz, ( 57
)
where

|E 30 | = lim |z|→∞ |E 3 |
is the asymptotic amplitude of the kink . Note that the availability of exact soliton solutions allows for the analytical calculation of the above integrals, hence of the condition (55), which is an extension of the well-known Vakhitov-Kolokov criterion. The equation (55) provides a sufficient stability condition [39]. Among all of the simultons belonging to the general family (52), the criterion (55) may only be applied within the specific constraint φ 1 = φ 2 . Indeed, in this case we found that the condition G = 0 leads, once again, to the previously found marginal stability condition

V = V cr = 2V 1 V 2 /(V 1 + V 2 )
. Moreover, a direct insight into the global stability properties of TWRIS for all possible values of their parameters may be gained with a powerful geometrical approach [42]. Indeed, TWRIS may be obtained as the solutions of the variational problem

δ(H + φ 1 I 12 -(φ 1 + φ 2 )I 23 -V J) = 0, ( 58
)
where δ is the Frechét derivative. In other words, TWRIS represent the extrema of the Hamiltonian (44), for a fixed value of the energies and momentum (here V, φ 1 , -(φ 1 + φ 2 ) represent Lagrange multipliers). Stable triplets are obtained whenever such extrema coincide with a global minimum of H. Clearly, if multiple solutions exist with the same I 12 , I 23 , J, the stable solution is obtained on the lower branch of H. In this framework, the condition (55) corresponds to solitons such that the normal vector to the three-dimensional surface H = H(I 12 , I 23 , J) lies in the space H = const. The above geometrical considerations permit the visualization of the stability boundaries when considering a projection of the hyper-surface H = H(I 12 , I 23 , J) on the plane (I 12 , H). For example, Fig. 5 displays the dependence of H upon I 12 for the case φ 1 = -φ 2 , where the criterion (55) cannot be applied, and in the case

φ 1 = φ 2 .
Here it is evident that the two branches of the Hamiltonian merge exactly at V = V cr : at this point, the normal to the H curve is also orthogonal to the vertical axis. Interestingly enough, Fig. 5 shows that the borderline TWRIS corresponds to a minimum of the bright pulses energy I 12 with respect to V . To summarize, we have shown that different numerical and analytical methods concur in predicting that the TWRIS stability is determined by the condition V > V cr .

Three-wave boomeron

The results of Fig. 4 illustrate that an unstable simulton with V < V cr decays into a stable simulton with V > V cr . This process is accompanied by the emission of an isolated bright pulse in wave E 3 . It is quite remarkable that the simulton decay and wave packet emission as it is numerically observed in Fig. 4 may be exactly reproduced in terms of an analytical higher-order soliton solution with varying speed, or "boomeron". Such solution was found by means of the techniques described in Ref. [15], and it can be expressed as

E 1 = 2pV 2 ∆ 2V 1 V 1 -V 2 e iq1z1 (H * + e -iθ -H * -e iθ ), (59)
E 2 = 2pV 1 ∆ 2V 2 V 1 -V 2 e iq2z2 (1 -Q)/(1 + Q)H + e i(β+θ) -(1 + Q)/(1 -Q)H -e -i(β+θ) , ( 60
)
E 3 = a V 1 V 2 e iq3z3 - ∆ 4p V 1 -V 2 V 1 V 2 E * 1 E * 2 , (61)
where

∆ = 1 + |H + | 2 1 + Q + |H -| 2 1 -Q -2 cos(β)Re(H + H * -e i(β+2θ) ) H ± (z, t) = e (-B±+iχ±)z e -2V 1 V 2 V 1 -V 2 (p-ik)t , ω = -2k V 1 V 2 V 1 -V 2 , χ ± = k V 1 + V 2 V 1 -V 2 ∓ 1 Q , B ± = p V 1 + V 2 V 1 -V 2 ∓ Q , tan(β) = k/(pQ), Q = 1 p 1 2 [ r + r 2 + 4k 2 p 2 ], r = p 2 -k 2 -a 2 , q n = q(V n+1 -V n+2 ), z n = z + V n t, n = 1, 2, 3 mod(3).
It should be pointed out that the expressions (59) and (60) of waves E 1 and E 2 could also be given in the form (17), that is in terms of the position ξ(t) and polarization vector α(t) of the boomeron by applying the transformations (41) and (42) to the expression (17). However, the equivalent expressions (59) and (60) are used here since they are more suitable to our analysis. Figure 6 displays the analytical boomeron solution corresponding to numerical dynamics of unstable simulton of Fig. 4. It is worth noting that the above solution depends upon seven real parameters V 1 , V 2 , p, k, q, θ, a > 0. From the definition of Q, it is apparent that these parameters must be chosen in such a way that if k = 0, then p 2 > a 2 . The analytical solution (59)-(61), while rather complicate at intermediate times, asymptotically consists of one or two coherent structures. In fact, let us consider first the decay process: if we assume p < 0, for negative large t (t → -∞) the boomeron is asymptotically composed of two bright pulses (E 1 , E 2 ) and a kink-like pulse (E 3 ) traveling with the locked velocity V i . If instead t is large and positive (t → +∞), the boomeron is composed of two bright pulses (E 1 , E 2 ) and a kink-like pulse (E 3 ) traveling at the locked velocity V f (V f > V i ), plus another pulse (E 3 ) that travels with the linear group velocity V 3 . The velocities V f and V i can be calculated from (59)-(61):

V i = 2V 1 V 2 V 1 + V 2 -Q(V 1 -V 2 ) (62) V f = 2V 1 V 2 V 1 + V 2 + Q(V 1 -V 2 ) . ( 63
)
The triplet traveling at very large |t| with the locked velocity V i (V f ) is in itself an exact solution of equations (40), namely it is the unstable (stable) TWRIS as presented in (48)-(50). Therefore, quite remarkably, the boomeron solution (59)-( 61) provides an exact description of the decay from unstable to stable simultons. Let us consider next the situation where a stable TWRIS collides with an isolated pulse in the wave E 3 , namely the excitation by absorption. Once again, this scattering process is exactly described by the boomeron solution ( 59)-( 61), and it leads to the excitation of an unstable TWRIS, which is induced by the absorption of the isolated wave E 3 . Indeed, whenever p > 0 and t is very large and negative, the boomeron (59)-( 61) is composed of a triplet consisting of two bright pulses (in waves E 1 , E 2 ) and a kink-like pulse (in wave E 3 ), all traveling with the same velocity V i , plus an isolated pulse in wave (E 3 ) that travels with the linear group velocity V 3 . The triplet and the isolated pulse collide and, as a result, the pulse in E 3 is completely absorbed by the triplet. Finally, for very large and positive t the boomeron consists of a single triplet formed by two bright pulses (in waves E 1 , E 2 ) and a kink-like pulse (in wave E 3 ), again traveling together with the velocity V f (V f < V i , see ( 62) and ( 63)). Note that the asymptotic boomeron triplets traveling with velocities V i and V f can be analytically mapped into the stable and unstable TWRIS as given in ( 48)-(50). In conclusion, the analytical solution (59)-( 61) with p > 0 provides the exact description of the excitation of an unstable TWRIS as a result of the inelastic collision between a stable TWRIS and an isolated wave packet. Figure 7(a) displays the analytical boomeron solution corresponding to the collision between a stable TWRIS and a pulse in wave E 3 . Whereas Fig. 7(b) shows the inelastic scattering of the TWRIS and the isolated pulse as numerically computed by integrating the equations (40) with the initial data at t = -5 equal to the solution of Fig. 7(a). As it can be seen in Fig. 7(b), the excited unstable TWRIS has a finite lifetime since it eventually decays into a stable or ground state TWRIS via the emission of another isolated wave. It is worth noting that both excitation and decay processes may be fully described by properly adjusting the parameters of Eqs.(59)-(61). The dynamics of the scattering between TWRIS and isolated waves is analogous to the interaction between a two-level atom and radiation. Indeed, transitions between excited and ground TWRIS states are induced by the absorption and spontaneous emission of an isolated pulse in the wave E 3 . Let us now briefly discuss the role of the various parameters in equations ( 59)-(61). Two of these parameters (i.e. the velocities V 1 and V 2 ) are fixed by the linear dispersive properties of the medium. We are thus left with five independent real parameters, namely p, k, q, a, θ (with the restrictions a > 0 and 0 ≤ θ < 2π). We would like to point out that, thanks to the above discussion, the specification of these parameters allows us to define all properties of both unstable and stable TWRIS. Indeed, these simultons result as asymptotic states of the boomeron expression (59)-(61) in the limit of |t| → ∞. The parameter p is associated with the rescaling of the wave amplitudes, and of the coordinates z and t. Whereas a measures the amplitude of the kink background in wave E 3 . The value of k is related to the soliton wave-number. The parameter q simply adds a phase shift which is linear in both z and t. Finally, θ fixes the shape of the stationary kink pulse E 3 . By adjusting the various degrees of freedom of the boomeron family of solutions ( 59)-(61), one may foresee the dynamical reshaping of the amplitude, phase, and velocity of the TWRIS. Additionally, one may fully describe the process of energy exchange among the three interacting waves.

Three-wave simulton collisions

In order to emphasize the novel and striking features of the scattering between TWRIS and isolated waves, let us briefly consider collisions between two different TWRIS. We denote with u n and v n (n = 1, 2, 3) the three components of the first and the second simulton, respectively, both of them centered at z = 0 for t = 0. Moreover, we will consider the following initial condition (at t = 0)

E 1 (z, 0) = u 1 (z - z 0 2 , 0) + v 1 (z + z 0 2 , 0)e iα , E 2 (z, 0) = u 2 (z - z 0 2 , 0) + v 2 (z + z 0 2 , 0)e -iα , E 3 (z, 0) = u 3 (z - z 0 2 , 0) × v 3 (z + z 0 2 , 0)/|u 3 (∞, 0)|,
where z 0 is the initial separation between the simultons, and α is a constant phase factor. Since equations (40) are completely integrable, one may intuitively expect (but not give for granted, in view of the previously discussed absorption and decay phenomena) that interactions between two initially well-separated TWRIS do not modify the shapes of triplets that emerge after the collision. Indeed, the numerical simulation of Fig. 8 shows that two TWRIS with different velocities penetrate and cross each other with no change of their shapes. The only effect of the interaction is a spatial shift and a phase shift, as it happens with ordinary bright TWRI solitons [43]. However, in a manner similar to cubic nonlinear Schrödinger solitons ( [44,[START_REF] Kodama | Soliton interaction in optical fibers[END_REF], [START_REF] Hasegawa | Solitons in Optical Communications[END_REF] and references therein), whenever the initial simulton separation is reduced, complex interaction phenomena may take place owing to the excitation of higher order soliton solutions. For example, Fig. 9(a) shows that two equal and in-phase (α = 0) TWRIS with the same velocity attract each other and periodically collapse. Whereas Fig. 9(b) shows that a repulsive force exists between two equal and out-of-phase (α = 0) solitons with the same velocity. In this case, two distinct TWRIS moving with different velocities emerge from the initial collision. Hence TWRI solitons may cross, attract or repel each other depending on their initial separation, velocity difference, and relative phase.

Conclusions and outlook

Nonlinear wave coupling may give rise to novel soliton structures such as the boomerons, whose dynamics is far richer than that of standard, uncoupled wave solitons. In this work we have shown that this is certainly the case for parametric resonant three-wave coupling, a model which is of wide applicability in a variety of physical contexts (plasma physics, nonlinear optics, acoustics, fluid dynamics, and solid state physics). In nonlinear optics, we envisage that propagation effects related to three-wave boomerons should be observable at realistic intensity levels (of the order of 10 M W/cm 2 ) when mixing picosecond light pulses in quadratic nonlinear materials. Indeed, it is already known that three-wave solitons may be exploited in order to enhance the efficiency of the frequency conversion of short optical pulses in the presence of group-velocity mismatch between the three waves [24,25,43]. Nevertheless, such process is hampered by the propagation instability of the up-converted pulse. This means that, for a given input pulse duration and intensity, a narrow range of tolerance results for the optimal nonlinear crystal length, thus limiting in practice the applicability of the soliton method. As we will discuss in separate reports, extended simultons and boomerons may provide a means to overcome this limitation: as it turns out, enhanced and stable parametric frequency conversion of short optical pulses may be achieved by mixing a short optical pulse with a quasi-continuous background pulse. As a further perspective, we believe that optical boomerons may open a new research direction in the growing field of slow and fast light, as they enable for the intensity-dependent and stable control of the speed of light pulses. Analytical solutions of the three-wave interaction equations (40) also point to the possibility of cavity-less short pulse train generation by means of the optical trappon effect (i.e., back and forth oscillations of self-trapped pulses inside a quadratic nonlinear medium, in the absence of any localized or distributed linearly reflecting mirror). We plan to investigate in the near future the possibility of generating boomeronic solutions in integrable stimulated scattering (i.e., Raman or Brillouin) three-wave equations, that is beyond the parametric mixing case of equations (40). As discussed in sections 1 and 2, we may anticipate that dynamically evolving solitons should exist and thus could be observed in a variety of nonlinear coupled mode evolution equations, including: parametric mixing in media with both quadratic and cubic nonlinearity, polarization evolution of short optical pulses in the presence of polarization mode dispersion and both instantaneous and relaxing nonlinearity (with a response time much longer than the pulse durations), etc. In conclusion, boomerons are definitely gaining the status of readily observable physical structures: their laboratory demonstration appears to be well at hand.

Fig. 1 .

 1 Fig. 1. Amplitude and phase of two simultons (48)-(50) at t = 0 with V1 = 2, V2 = 1, a = 1, k = 0.5, q = 0.5; in the upper figures p = 1, in the lower figures p = -1. Thin solid curve: E1; dashed curve: E2; thick solid curve: E3.

Fig. 2 .

 2 Fig. 2. Mapping of parameters (V, φ1), in the space φ1 = -φ2, into parameters (p, k) (q = 0, a = 1), for V1 = 2, V2 = 1. Shaded areas represent forbidden parameter regions.

Fig. 3 .

 3 Fig. 3. (a) phase constant φ1 versus velocity V for different energies I12. (b) energies I1, I2 versus velocity V for representative values of the phase constant φ1. (c) energies I12 versus velocity V for representative values of the phase constant φ1. (d) transverse momentum J versus velocity V . In all cases φ2 = -φ1, V1 = 2, V2 = 1.

3 Fig. 4 .

 34 Fig.4. Propagation of stable and unstable simultons. The common parameters are V1 = 2, V2 = 1, a = 1, k = 0.5, q = 1. For the stable simulton (upper figures) p = 1 (V = 1.68 > Vcr, φ1 = 0.5, φ2 = 2.2), whereas for the unstable simulton (lower figures) p = -1 (V = 1.1 < Vcr, φ1 = 1.43, φ2 = -0.33).

2 Fig. 5 .

 25 Fig. 5. (a) Hamiltonian H versus energy I12 and (b) I12 versus velocity V . Characteristic velocities are V1 = 2, V2 = 1 and φ2 = -0.2. Dashed (solid) curves correspond to unstable (stable) solitons.

Fig. 6 .

 6 Fig. 6. Analytical boomeron solution describing the unstable TWRIS dynamics of Fig. 4. Parameters are V1 = 2, V2 = 1, V3 = 0, p = -1, a = 1, k = 0.5, q = 1, θ = π/2 (The triplet velocities are Vi = 1.1 and V f = 1.68 (Vcr ≈ 1.3)).

Fig. 7 .

 7 Fig. 7. a) Analytical boomeron solution describing the collision of a stable TWRIS with a single pulse in wave E3. Parameters are V1 = 2, V2 = 1, V3 = 0, p = 1, a = 1, k = 0.5, q = 1, θ = π/6. The triplet velocities are Vi = 1.68 and V f = 1.1 (Vcr ≈ 1.3). b) Numerical excitation-decay process.

Fig. 8 .

 8 Fig. 8. Collision of two stable TWRIS with different velocities. Fast simulton V = 1.9, slow simulton V = 1.7. Simulation is performed in reference frame moving at velocity V ref = 1.8.

Fig. 9 .

 9 Fig. 9. a) Collision of two equal and in-phase stable TWRIS with the same velocity V = 1.8; b) collision of two equal and π out-of-phase stable TWRIS with the same velocity. Simulations are performed in reference frame moving at velocity V ref = 1.8.


Ph.D S Randall Thomas

Mathematical models for kidney function focusing on clinical interest (and in French: Modèles mathématiques du rein dans l'intérêt du clinicien)

Keywords: VPH, Mathematical models, kidney anatomy, renal physiology, electrolyte metabolism, renal epithelial transport ------------------ VPH, Modèles mathématiques, anatomie rénale, physiologie rénale, métabolisme des électrolytes, transport épithélial rénal

<= 250 words) We give an overview of mathematical models of renal physiology and anatomy with the clinician in mind. Beyond the past focus on issues of local transport mechanisms along the nephron and the urine concentrating mechanism, recent models have brought insight into difficult problems such as renal ischemia (oxygen and CO2 diffusion in the medulla) or calcium and potassium homeostasis. They have also provided revealing 3D reconstructions of the full trajectories of families of nephrons and collecting ducts through cortex and medulla. The recent appearance of sophisticated whole-kidney models representing nephrons and their associated renal vasculature promises more realistic simulation of renal pathologies and pharmacological treatments in the foreseeable future.

Introduction

In this review, we summarize mathematical modelling studies of kidney physiology and anatomy that have relevance for the clinician. We make no pretense of giving the detailed contributions or modeling techniques of the mentioned modeling studies, but rather point to the relevant studies for those interested in particular issues. For more general reviews of renal modelling, we refer the reader to the following papers:

• on the urine concentrating mechanism, see Stephenson [START_REF] Stephenson | Urinary concentration and dilution: models[END_REF] for a thorough, technically detailed history starting with Hargitay & Kuhn [START_REF] Hargitay | Das Multipikationsprinzip als Grundlage der Harnkonzentrierung in der Niere[END_REF][START_REF] Hargitay | The multiplication principle as the basis for concentrating urine in the kidney[END_REF], and see Pannabecker and colleagues for more recent developments that take advantage of imaging studies and 3D reconstruction of medullary nephron-vascular anatomy [START_REF] Dantzler | Urine-concentrating mechanism in the inner medulla: function of the thin limbs of the loops of Henle[END_REF][START_REF] Neal | Novel hemodynamic structures in the human glomerulus[END_REF][START_REF] Stephenson | Concentration of urine in a central core model of the renal counterflow system[END_REF]; • for reviews and models of tubular function and tubule-specific transporters, see

Weinstein [START_REF] Weinstein | Mathematical models of tubular transport[END_REF][START_REF] Weinstein | Insights from mathematical modeling of renal tubular function[END_REF][START_REF] Weinstein | Mathematical models of renal fluid and electrolyte transport: acknowledging our uncertainty[END_REF][START_REF] Weinstein | A Mathematical Model of the Rat Nephron: Glucose Transport[END_REF][START_REF] Weinstein | A Mathematical Model of Rat Proximal Tubule and Loop of Henle[END_REF][START_REF] Weinstein | Systems biology of the cortical collecting duct[END_REF][START_REF] Weinstein | A mathematical model of the rat kidney: K(+)-induced natriuresis[END_REF]; and • for more general reviews of renal modelling at several scales, see Edwards [START_REF] Edwards | Modeling transport in the kidney: investigating function and dysfunction[END_REF],

Thomas et al. [START_REF] Thomas | Kidney modelling: status and perspectives[END_REF], and Thomas [START_REF] Thomas | Kidney Modeling and Systems Physiology[END_REF]. We also point out that experimenters and clinicians are well advised to read the modelling literature -the models are often the best sort of critical review of the experimental literature, since they must establish a coherent set of model parameters based on experimental measurements made under a variety of conditions using a number of experimental techniques and often in several different species.

Models of whole-kidney function

Among the earliest models of whole-kidney function are the treatments by Guyton and colleagues. Developed as part of the classic Guyton model of cardiovascular function [START_REF] Guyton | Systems analysis of arterial pressure regulation and hypertension[END_REF][START_REF] Guyton | Circulation: overall regulation[END_REF], these were not stand-alone kidney models and they gave a rather simple input-output description of kidney function with scant internal detail, but they were developed with the clinic in mind by Guyton and his colleagues and established their conviction that the kidney plays a central role in blood pressure regulation. This work was of central importance for research into the mechanisms of blood pressure regulation, though a recent paper [START_REF] Kurtz | Testing Computer Models Predicting Human Responses to a High-Salt Diet[END_REF] has cast doubt on the validity of the predictions of both the classic Guyton model and its more recent descendants [START_REF] Abram | Quantitative Circulatory Physiology: an integrative mathematical model of human physiology for medical education[END_REF][START_REF] Hester | HumMod: A Modeling Environment for the Simulation of Integrative Human Physiology[END_REF] under conditions of salt loading. Hallow and colleagues [START_REF] Hallow | Primary proximal tubule hyperreabsorption and impaired tubular transport counterregulation determine glomerular hyperfiltration in diabetes: a modeling analysis[END_REF][START_REF] Hallow | Evaluation of renal and cardiovascular protection mechanisms of SGLT2 inhibitors: model-based analysis of clinical data[END_REF][START_REF] Hallow | A model-based approach to investigating the pathophysiological mechanisms of hypertension and response to antihypertensive therapies: Extending the Guyton model[END_REF] have updated and improved on the Guyton "platform", taking advantage of Karaaslan's improved version [START_REF] Karaaslan | A mathematical model of long-term renal sympathetic nerve activity inhibition during an increase in sodium intake[END_REF][START_REF] Karaaslan | Long-term mathematical model involving renal sympathetic nerve activity, arterial pressure, and sodium excretion[END_REF] of the Guyton kidney. Moss et al. [START_REF] Moss | Virtual patients and sensitivity analysis of the Guyton model of blood pressure regulation: towards individualized models of whole-body physiology[END_REF] presented a scenario for future individualized, "virtual patient" models of whole-body physiology based on an extension of the classic Guyton model.

At the other extreme of detail are models treating solute and water transport along all the nephrons in the whole kidney or in the whole medulla. Following Weinstein [START_REF] Weinstein | A mathematical model of the rat kidney: K(+)-induced natriuresis[END_REF], one can distinguish nephron models from kidney models -in a nephron model, peritubular or interstitial conditions are specified (i.e., they are boundary conditions, not variables), whereas in a kidney model, transport by tubules and blood vessels determine the peritubular environment. Of kidney models, Weinstein's recent model [START_REF] Weinstein | A mathematical model of the rat kidney: K(+)-induced natriuresis[END_REF] is the most detailed, since it is the first to include the vasculature (vasa recta) in a model that treats epithelial transport of a rich set of ions and buffer systems along each nephron segment in the whole population of short and long nephrons in both cortex and medulla. This is clearly the most complete whole-kidney model to date, and it is replete with cogent reflection on many issues of clinical relevance. The whole-kidney model by Moss & Thomas [41] also includes both nephrons and vasculature but treats tubular epithelia as a simple barrier; their A c c e p t e d M a n u s c r i p t model was the first to explicitly couple glomerulovascular and medullary dynamics and aims to serve as a platform for simulation of pathophysiological conditions and pharmacological interventions. The Moss & Thomas [START_REF] Moss | Hormonal regulation of salt and water excretion: a mathematical model of whole kidney function and pressure natriuresis[END_REF] model was also conceived to be embedded in whole body physiological models for exploration of the relationship between the extrarenal environment and renal excretory function in physiological and pathophysiological contexts.

The three previous families of kidney models treated nephron-vascular exchanges of only salt and urea and were focused on the still open question of the origin of the corticomedullary gradient of osmolality and the urine concentrating mechanism (see below). These three families of models were those by Stephenson and colleagues [START_REF] Stephenson | Urinary concentration and dilution: models[END_REF][START_REF] Stephenson | Electrolyte, urea, and water transport in a twonephron central core model of the renal medulla[END_REF], by Layton and Layton [START_REF] Layton | A region-based mathematical model of the urine concentrating mechanism in the rat outer medulla. I. Formulation and base-case results[END_REF][START_REF] Layton | Countercurrent Multiplication May Not Explain the Axial Osmolality Gradient in the Outer Medulla of the Rat Kidney[END_REF], and that of Wexler and colleagues [START_REF] Wexler | Three-dimensional anatomy and renal concentrating mechanism. I. Modeling results[END_REF], advanced by Thomas [START_REF] Hervy | Inner medullary lactate production and urine-concentrating mechanism: a flat medullary model[END_REF][START_REF] Thomas | Cycles and separations in a model of the renal medulla[END_REF][START_REF] Thomas | Inner medullary external osmotic driving force in a 3-D model of the renal concentrating mechanism[END_REF][START_REF] Wang | Outer medullary anatomy and the urine concentrating mechanism[END_REF].

Future modeling work will profit from new insights from experimental manipulations, such as recent work on "anomalous" tubuloglomerular feedback (TGF) in rats fed a high sodium diet, such that stimulation of TGF increases rather than decreases single nephron GFR, thereby favoring excretion of a high sodium load rather than stabilization of GFR [START_REF] Braam | Overruled": the kidneys' judgment of sodium balance versus stabilization of renal function[END_REF][START_REF] Thomson | Nitric oxide mediates anomalous tubuloglomerular feedback in rats fed high-NaCl diet after subtotal nephrectomy[END_REF].

Renal epithelia and specificities of tubular function

Transport along each segment of the nephron is specific to the functions local to the given segment, and each segment has specific membrane transporters and epithelial permeability properties. Unraveling the specific roles of these transporters and the pathologies arising from transporter defects has profited greatly from give and take between experimental and modeling studies. These epithelial and tubular models [59-62, 65,66, 2,9-11,34] have provided much clinical insight and are generally motivated not by a desire to simply simulate published experimental results but rather to resolve conflicting interpretations of experimental data or address open questions and clinical conundrums.

Potassium handling

Plasma potassium concentration ([K + ]) must be strictly regulated, since even slight elevation above normal levels incurs serious cardiac and neurological disorders. The kidney is the principal regulator of K + homeostasis [START_REF] Meneton | Sodium and potassium handling by the aldosterone-sensitive distal nephron: the pivotal role of the distal and connecting tubule[END_REF], and the site of this regulation is the distal nephron, where K + secretion is managed. Several modeling studies [START_REF] Layton | Renal potassium handling in rats with subtotal nephrectomy: modeling and analysis[END_REF][START_REF] Weinstein | A mathematical model of rat distal convoluted tubule (II): Potassium secretion along the connecting segment[END_REF][START_REF] Weinstein | Potassium excretion during antinatriuresis: perspective from a distal nephron model[END_REF][START_REF] Weinstein | A mathematical model of the rat kidney: K(+)-induced natriuresis[END_REF][START_REF] Yang | Regulation of renal Na transporters in response to dietary K[END_REF] have addressed questions of renal potassium handling. In a model of epithelial transport along the full distal nephron, Weinstein considers the question of independent regulation of Na + and K + excretion, focusing on familial hyperkalemic hypertension (FHH), concluding that there must be two distinct distal nephron signals. The recent whole-kidney model of Weinstein [START_REF] Weinstein | A mathematical model of the rat kidney: K(+)-induced natriuresis[END_REF] (see above) provides the most complete context for quantitative speculation about the renal responses to perturbations of (among other things) electrolyte balance. In this paper, one issue explored in some depth is the cascade of effects on transport of sodium and potassium in proximal tubule, ascending Henle's limb, and distal tubule resulting from mild hyperkalemia (plasma K + was raised from 5.0 to 5.5 mM) -the model identifies a 3.3% decrease of proximal tubule Na + reabsorption as the source of the predicted 2.3-fold increase in Na + excretion (K + excretion increases by 1.3-fold). Layton, Edwards, & Vallon [START_REF] Layton | Renal potassium handling in rats with subtotal nephrectomy: modeling and analysis[END_REF], in their whole-nephron model (but with no explicit vasculature), directly addressed renal adaptations to increases of K + load under normal conditions or following reduced renal function. They quantitatively evaluated the limits to excretion of an acute or chronic potassium load in simulations of nephrectomized rats.

A c c e p t e d M a n u s c r i p t

Calcium handling

As reviewed by Friedman [START_REF] Friedman | Codependence of renal calcium and sodium transport[END_REF][START_REF] Friedman | Mechanisms of renal calcium transport[END_REF] and Moor & Bonny [START_REF] Moor | Ways of calcium reabsorption in the kidney[END_REF], overall calcium reabsorption in the kidney proceeds in parallel with sodium reabsorption, especially in the proximal tubule and loop of Henle (with different cellular mechanisms). However, hormonal and pharmacological regulation of renal calcium balance is situated in the distal nephron, where, under the influence of PTH and calcitonin and other factors, sodium and calcium transport can be dissociated and are inversely related. The channels and transporters involved in the various cell types along the nephron have (largely) been worked out, including identification of the mutations responsible for inherited disorders of distal calcium absorption such as Bartter's and Gitelman's syndromes. Nonetheless, questions remain, and mathematical models are needed to understand quantitatively the details of the kidney's handling of calcium transport as a function of hormonal signals, blood calcium and phosphate levels, and the various factors that perturb these.

To this purpose, recent mathematical models of calcium transport along the proximal tubule [START_REF] Edwards | A model of calcium transport and regulation in the proximal tubule[END_REF] and along the late nephron [START_REF] Bonny | Calcium reabsorption in the distal tubule: regulation by sodium, pH, and flow[END_REF][START_REF] Edwards | Regulation of calcium reabsorption along the rat nephron: a modeling study[END_REF] (thick ascending limb, distal tubule, cortical and medullary collecting ducts) have contributed quantitative and qualitative substance to these issues. In particular, building on a model of distal tubule calcium transport [START_REF] Bonny | Calcium reabsorption in the distal tubule: regulation by sodium, pH, and flow[END_REF] and an integrated epithelial model of ion transport from the thick ascending limb through the inner medullary collecting duct [START_REF] Weinstein | A mathematical model of rat ascending Henle limb. III. Tubular function[END_REF], Edwards [START_REF] Edwards | Regulation of calcium reabsorption along the rat nephron: a modeling study[END_REF] was the first to include calcium transport in a detailed model of epithelial transport from the medullary thick ascending limb through the inner medullary collecting duct. Her model focused especially on the role of calcium-sensing mechanisms (e.g., CaSR, the calcium-sensing receptor). For calcium transport in the proximal convoluted tubule (PCT), a recent model by Edwards & Bonny [START_REF] Edwards | A model of calcium transport and regulation in the proximal tubule[END_REF] confirmed quantitatively that the generally accepted discursive model of transcellular and paracellular calcium transport is consistent with available experimental data (principally from micropuncture studies in rodents), and provided detailed predictions of the expected impact of PTH, 1,25(OH)2D3, and calcitonin or perturbations of cell membrane transporters (NHE3 for protons and sodium, SGLT2 for sodium-glucose) or tight-junction permeabilities.

As pointed out by Weinstein [START_REF] Weinstein | By the way, proximal tubule calcium transport[END_REF], the models of Edwards and Bonny [START_REF] Bonny | Calcium reabsorption in the distal tubule: regulation by sodium, pH, and flow[END_REF][START_REF] Edwards | Regulation of calcium reabsorption along the rat nephron: a modeling study[END_REF][START_REF] Edwards | A model of calcium transport and regulation in the proximal tubule[END_REF] provide the key components for an eventual overall treatment of renal calcium and phosphate handling. The next step would be to combine them into an integrated model (more complete and detailed than the one by Granjon et al. [START_REF] Granjon | A model of calcium homeostasis in the rat[END_REF]) that would include PCT, short and long loops of Henle, and distal nephrons of both superficial and juxtamedullary nephrons to enable simulation of acute changes in kidney function, blood composition, or diuretic administration.

Oxygen delivery/consumption in cortex and medulla…

The kidney should be the organ that is least susceptible to hypoperfusion and ischemia since it has by far the highest rate of blood flow per gram of tissue (420 ml/min/100g in kidney, but 54 ml/min/100g in brain and 84 ml/min/100g in heart) and its ratio of O2 consumption/ O2 delivery is only about 8%, which is lower than that for any other organ (cf. 34% for brain and skeletal muscle and 65% for heart) (values in humans, from Table 1 of Brezis et al. [START_REF] Brezis | Renal ischemia: a new perspective[END_REF]). Nonetheless, the kidney suffers the most severe hypoperfusion damage in patients with acute renal failure (ARF). The explanation for this seeming paradox is two-fold. First, oxygenation and blood flow are heterogeneous within the kidney due to the essentially separate irrigation of cortex, outer medulla, and inner medulla and also to the special hairpin-turn configuration of nephrons and vessels (and consequent recycling and counter-A c c e p t e d M a n u s c r i p t current exchange of solutes and water). And secondly, transepithelial solute transport consumes a high amount of energy, particularly in the mTAL (medullary thick ascending limb), situated in the relatively hypoxic region of the outer medulla outside the vascular bundles -for outer medulla, the ratio of O2 consumption/ O2 delivery is at least 79% (rat kidney, Table 1 of Brezis et al. [START_REF] Brezis | Renal ischemia: a new perspective[END_REF]). As Brezis et al. [START_REF] Brezis | Renal ischemia: a new perspective[END_REF] also pointed out, decreased urine concentrating ability (which depends on mTAL active transport) is the most consistent defect in ARF (acute renal failure), is a very early sign of impending ARF, and is exquisitely sensitive to brief periods of experimental renal hypoxia induced by cyanide or mild renal ischemia even in the absence of ARF.

Thus, given that the special problems of renal oxygenation are intimately tied up with the complex structure of the kidney, and despite the technical difficulties (high diffusion rates, complex buffering relations, limited experimental measurements of certain important parameters, etc.), a number of mathematical models of renal oxygen transport and metabolism have been realized, ranging in complexity from a simple but clinically motivated Matlab model [START_REF] Burke | Factors maintaining a pH gradient within the kidney: role of the vasculature architecture[END_REF] to sophisticated medullary [START_REF] Chen | A mathematical model of O2 transport in the rat outer medulla. I. Model formulation and baseline results[END_REF][START_REF] Fry | Impact of renal medullary three-dimensional architecture on oxygen transport[END_REF][START_REF] Zhang | Oxygen transport across vasa recta in the renal medulla[END_REF] or even whole kidney 3D models [START_REF] Layton | A computational model for simulating solute transport and oxygen consumption along the nephrons[END_REF][START_REF] Lee | A model of oxygen transport in the rat renal medulla[END_REF] of oxygen delivery, distribution, and consumption under normal and disturbed conditions. These modeling studies and related experimental studies firmly establish the delicate balance maintained by the kidney between, on one hand, minimizing medullary blood flow in order to favor increased concentrating ability, and on the other hand, increasing medullary blood flow to deliver sufficient oxygen for the metabolically expensive transport processes of the mTAL and the medullary CD. These models also explore various questions relating to hypothesized tradeoffs between aerobic and anaerobic modes of transport within the low PO2 environments of outer and inner medulla.

Urine concentrating mechanism

As pointed out above, abnormal dilution or concentrating ability of the urine is often an early sign of renal pathologies, and while the main principles of the concentrating mechanism are well understood for the outer medulla (see, e.g., [START_REF] Layton | Urine concentrating mechanism: impact of vascular and tubular architecture and a proposed descending limb urea-Na+ cotransporter[END_REF][START_REF] Layton | Permeability criteria for effective function of passive countercurrent multiplier[END_REF][START_REF] Thomas | Kidney Modeling and Systems Physiology[END_REF]), the origin of the steep osmotic gradient in the inner medulla remains one of the central mysteries of renal physiology. For example, it remains unknown how the kidney is able to concentrate urine to the degree that is evident in antidiuresis, and so this cannot be modeled mechanistically without hypothesizing an extra driving force such as metabolic osmole production in the medullary interstitium (perhaps lactate production by anaerobic glycolysis [START_REF] Hervy | Inner medullary lactate production and urine-concentrating mechanism: a flat medullary model[END_REF][START_REF] Jen | Externally driven countercurrent multiplication in a mathematical model of the urinary concentrating mechanism of the renal inner medulla[END_REF][START_REF] Thomas | Inner medullary lactate production and accumulation: A vasa recta model[END_REF][START_REF] Thomas | Inner medullary external osmotic driving force in a 3-D model of the renal concentrating mechanism[END_REF]) or some other mechanism [START_REF] Knepper | Concentration of solutes in the renal inner medulla: interstitial hyaluronan as a mechano-osmotic transducer[END_REF][START_REF] Layton | Urine concentrating mechanism: impact of vascular and tubular architecture and a proposed descending limb urea-Na+ cotransporter[END_REF]. This ongoing story is too involved to relate here, but has been described in a number of reviews [START_REF] Dantzler | Urine-concentrating mechanism in the inner medulla: function of the thin limbs of the loops of Henle[END_REF][START_REF] Edwards | Modeling transport in the kidney: investigating function and dysfunction[END_REF][START_REF] Pannabecker | Role of three-dimensional architecture in the urine concentrating mechanism of the rat renal inner medulla[END_REF][START_REF] Stephenson | Urinary concentration and dilution: models[END_REF][START_REF] Thomas | Inner medullary lactate production and accumulation: A vasa recta model[END_REF][START_REF] Thomas | Kidney Modeling and Systems Physiology[END_REF], and most recently by Nawata & Pannabecker [START_REF] Nawata | Mammalian urine concentration: a review of renal medullary architecture and membrane transporters[END_REF], who relate in particular the collaborative contributions of anatomical and functional modeling.

Models of renal anatomy

Most mathematical models of renal anatomy have been reconstructions from serial sections, focused on the medulla or on the glomerulus. Nordsletten [START_REF] Nordsletten | Structural morphology of renal vasculature[END_REF] gave a 3D reconstruction (see Figure [START_REF] Abram | Quantitative Circulatory Physiology: an integrative mathematical model of human physiology for medical education[END_REF]) of arterial and venous vessels of the whole rat kidney down to 30 µm in diameter and made the 3D dataset available for download for reconstruction and morphologically based blood flow analysis. Pannabecker and colleagues have published a number of studies reconstructing nephrons and collecting ducts in the rat medulla (for example, Figure [START_REF] Bonny | Calcium reabsorption in the distal tubule: regulation by sodium, pH, and flow[END_REF]); these are reviewed by Newata & Pannabecker [START_REF] Nawata | Mammalian urine concentration: a review of renal medullary architecture and membrane transporters[END_REF] in the context of recent advances towards understanding the urine concentrating mechanism. Zhai, Christensen, and colleagues published wonderfully revealing reconstructions of the paths of multiple families of nephrons and vessels in mouse [START_REF] Letts | Nephron morphometry in mice and rats using tomographic microscopy[END_REF][START_REF] Zhai | Three-Dimensional Reconstruction of the Mouse Nephron 10.1681/ASN.2005080796[END_REF][START_REF] Zhai | Digital Three-Dimensional Reconstruction and Ultrastructure of the Mouse Proximal Tubule[END_REF] and rat [START_REF] Christensen | Threedimensional reconstruction of the rat nephron[END_REF][START_REF] Letts | Nephron morphometry in mice and rats using tomographic microscopy[END_REF] A c c e p t e d M a n u s c r i p t kidneys, with great relevance for issues of nephron-vascular exchange and the urine concentrating mechanism. These have provided great insight into the intimate details of the relative placement of long and short nephrons and vasculature (vasa recta, in particular), relevant especially for questions of salt and urea recycling in studies of the urine concentrating mechanism.

A c c e p t e d M a n u s c r i p t

Neal et al. [START_REF] Neal | Novel hemodynamic structures in the human glomerulus[END_REF][START_REF] Neal | Three-dimensional reconstruction of glomeruli by electron microscopy reveals a distinct restrictive urinary subpodocyte space[END_REF][START_REF] Neal | Glomerular filtration into the subpodocyte space is highly restricted under physiological perfusion conditions[END_REF] published 3D reconstructions of glomeruli, not only in the rat but also in human kidneys. They were the first to show [START_REF] Neal | Three-dimensional reconstruction of glomeruli by electron microscopy reveals a distinct restrictive urinary subpodocyte space[END_REF] that Bowman's Space is actually three functionally distinct spaces. They describe [START_REF] Neal | Novel hemodynamic structures in the human glomerulus[END_REF] the functional implications of the difference in scale between rodent and human kidneys. Their most recent studies have clarified the roles of these spaces in glomerular filtration and point to possible correlations with glomerular pathologies.

Other models of practical clinical use

Other "models", which may not qualify as models of the kidney, can be useful didactic tools for clinical procedures. An example is Walther et al. [START_REF] Walther | Downloadable computer models for renal replacement therapy[END_REF], which provides two spreadsheet calculators for prediction of solute clearances and solute concentrations during renal replacement therapy -these are useful for teaching but not for clinical application, since they have not been clinically tested or FDA approved. The first calculates clearance values from inputs including blood flow, hematocrit, ultrafiltration rate, dialysate flow rate, as well as the reflection coefficient and the mass transfer area coefficient for the solute of interest, and the free fraction if the solute is protein bound. The second program calculates plasma solute concentrations from the clearance values obtained by the first program and from additional input values including the number of treatments per week, the duration of the treatments, and the solute's production rate and volumes of distribution. The program calculates the time-averaged solute concentration and provides a graphic display of the solute concentration profile through a week-long interval.

Conclusions

Given its unique triple role of maintaining water balance, salt balance, and nitrogen balance, the mammalian kidney fulfills such a multitude of functions that it would be vain to expect any single mathematical model to capture them all. As shown here, many sorts of models have been published, from simple black-box treatments of overall renal input-output relations to extremely sophisticated models of water and electrolyte transport across the epithelia of all nephrons and vessels. Certainly, the great majority of these models treat the rat kidney, since it is for rats and mice that there is the most complete set of experimental measurements. Nonetheless, the motivation for most of these studies is to further the understanding of human kidney function, and in general they carefully relate the model results to issues of human physiology and pathophysiology. One can therefore wholeheartedly encourage the clinician to delve into the modeling literature despite its necessarily technical nature -each mathematical modeling paper is a carefully reasoned review of the relevant experimental literature (two outstanding examples are Moss & Thomas [START_REF] Moss | Hormonal regulation of salt and water excretion: a mathematical model of whole kidney function and pressure natriuresis[END_REF] and Weinstein [START_REF] Weinstein | A mathematical model of the rat kidney: K(+)-induced natriuresis[END_REF]).

Fig. 1 .

 1 Fig. 1. Side view of rat renal vasculature with vessel thickness proportional to a quarter the true radial value. Color of the arterial tree indicates vessel order stability, white giving {rho}s = 0 and red giving {rho}s = 1. The venous tree exhibits the same information on a fading blue scale. (from reference [46]: Nordsletten et al. 2006).

Fig. 2 .

 2 Fig. 2. Three-dimensional reconstruction of all papillary CDs (AQP2, blue), ATLs (ClC-K1, green), and DTLs ({alpha} B-crystallin, yellow) from a single kidney. Papillary surface epithelium is shown in gray. A: lateral view. B: axial view from papilla tip. Scale bars = 100 µm. See Supplemental Video 2 for additional details.( from reference [47], Pannabecker and Dantzler 2007.)

Figure 3

 3 Figure 3 near here Fig. 3. All SLN, LLN, and CD traced within one kidney. The red lines represent the renal capsule and the papilla. The white dots represent glomeruli (not actual size). Tubule segments are represented in different colors: PT, blue; thin limbs (TL), green; TAL, red; DCT, purple; connecting tubules (CNT), orange; and CD, brown.(B) Level of the OSOM (approximately 1570 µm from the renal surface). (C) Level of the initial part of the ISOM (approximately 1925 µm). (D) Level of the middle to lower ISOM (approximately 2585 µm). White circles highlight the same vascular bundles. . Scale = 50 µm. (from[START_REF] Zhai | Three-Dimensional Reconstruction of the Mouse Nephron 10.1681/ASN.2005080796[END_REF], reference[START_REF] Zhai | Three-Dimensional Reconstruction of the Mouse Nephron 10.1681/ASN.2005080796[END_REF])

A c c e

  p t e d M a n u s c r i p t A c c e p t e d M a n u s c r i p t

Acknowledgements

The author is grateful to the CNRS for continuing support during his emeritus period.

Conflict of interest

No potential conflict of interest relevant to this article was reported.