UNIVERSITÀ DEGLI STUDI LA
SAPIENZA
Facoltà di Scienze Matematiche, Fisiche e Naturali
Corso di Laurea Specialistica in Fisica
Teoria di Whitham per lo studio di
shock dispersivi in sistemi non
integrabili
Relatori:
Prof. Giancarlo Ruocco
Prof. Andrea Fratalocchi
Prof. Stefano Trillo
Anno Accademico 2010-2011
Candidato:
Matteo Crosta
1
...Ai miei gentori
Abstract
This thesis investigated the dynamics of dispersive shocks in nonintegrable
systems of physical interest. The study is carried out by employing a fully analytic framework given by the Whitham theory of modulation. This approach,
originally developed by Whitham to address the shocks dynamics of integrable
equations, has been recently extended to nonintegrable systems thus providing
a powerful tool to analytically study dispersive shocks in a broad range of situations. In the first part of the thesis, I presented a deep introduction to the
Whitham theory, discussing both integrable and nonintegrable equations. I then
applied this theory to study the dispersive shock dynamics of the Cubic-Quintic
Nonlinear Schrödinger Equation (CQNLS), an equation of utmost importance in
both classical physics and quantum mechanics. More specifically, I studied the
development of a dispersive shock from a sharp, step-like initial condition, which
represents a bright pulse in optics or a localized matter-wave in Bose-Einstein
Condensate (BEC). I begin the analysis by deriving a complete analytic set of
periodic solutions of the CQNLS, and discussed their stability and bistability
properties. The results of this work have fostered a first publication [1]. I then
focused the analysis on the dispersive shock, and through the Whitham theory
I carried out a fully analytic study of the phenomenon [2].
This thesis led to the following publications:
[1] M. Crosta, A. Fratalocchi and S. Trillo, ”Bistability and instability of
dark-antidark solitons in the cubic-quintic nonlinear Schrödinger equation”,
Phys. Rev. A (in review).
[2] M. Crosta, S. Trillo and A. Fratalocchi, ”Dispersive shocks in BoseEinstein Condensates: dynamics of long range interactions”, Phys. Rev. Lett.
(in preparation).
Indice
1 Introduzione
2
2 Teoria di Whitham
2.1 Formazione di shock . . . . . . . . . . . . . . . . . .
2.2 Teoria di Whitham: formulazione generale . . . . . .
2.3 Teoria di Whitham nella KdV . . . . . . . . . . . . .
2.3.1 Onde Cnoidali . . . . . . . . . . . . . . . . .
2.3.2 Equazioni di modulazione . . . . . . . . . . .
2.3.3 Diagonalizzazione delle leggi di conservazione
2.4 Costruzione di Gurevich-Pitaevskii . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4
4
6
8
8
12
16
20
3 Teoria di Whitham in sistemi non integrabili
3.1 Formulazione non integrabile per la KdV . . . . . . . .
3.1.1 Trailing Edge . . . . . . . . . . . . . . . . . . .
3.1.2 Leading Edge . . . . . . . . . . . . . . . . . . .
3.2 Formulazione non integrabile per sistemi bidirezionali
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
25
25
27
30
33
4 Nonlinear Schrödinger quintica
36
4.1 Teoria analitica dello shock generato da una discontinuità iniziale 37
4.1.1 Limite non dispersivo . . . . . . . . . . . . . . . . . . . . 38
4.2 Solitoni Dark e Antidark . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.1 Bistabilità e instabilità . . . . . . . . . . . . . . . . . . . . 48
4.2.2 Parametrizzazione della soluzione di traslazione negli shock
dispersivi . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3 Shock dispersivo . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.3.1 Leading Edge . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3.2 Trailing Edge . . . . . . . . . . . . . . . . . . . . . . . . . 57
5 Conclusioni
61
A Invarianti di Riemann
63
B Calcolo di onde cnoidali per la NLS Quintica
67
1
Capitolo 1
Introduzione
Gli shock sono fenomeni fisici legati a rapidi e quasi discontinui cambiamenti
delle proprietà di un mezzo, come pressione e densità di un fluido, dove si propaga un onda detta ”shock wave”. Questi fenomeni sono legati esclusivamente
ad effetti non lineari. In questi sistemi, dopo un tempo finito, le condizioni
iniziali evolvono in singolarità formando una regione in cui la soluzione è ”non
fisica”. Tale comportamento si regolarizza attraverso l’aggiunta di piccoli termini di diffusione o dispersione che diventano rilevanti per tempi confrontabili
con il tempo di formazione della singolarità. A seconda di che tipo di regolarizzazione si aggiunge nel modello, la fisica varia drasticamente. Nei mezzi
viscosi, ad esempio gas comprimibili, si parla di viscous shock wave (VSW),
caratterizzate da uno smorzamento del fronte d’onda che evita la formazione
della regione singolare. Nei mezzi dispersivi, invece, si osserva la formazione di
una intensa regione oscillatoria che tipicamente termina in un solitone a seguito
della compensazione tra la dispersione e la non linearità. In questo caso si parla
di dispersive shock wave (DSW).
In questa tesi verrà focalizzata l’attenzione sui modelli dispersivi e i fenomeni di shock ad essi associati. Nella regione oscillatoria dei DSW si distinguono
tipicamente due diverse scale di evoluzione: la scala veloce di oscillazione del
treno d’onda e la scala lenta di modulazione del profilo oscillatorio. Nel 1965
Whitham [1] formulò un metodo di mediazione per equazioni differenziali non
lineari a derivate parziali (PDEs). Il metodo proposto è una generalizzazione
del teorema della media per equazioni differenziali che ammettono una soluzione
di propagazione del tipo f (x − vt). L’idea generale è che mediando le leggi di
conservazione rispetto alla fase veloce θ = x − vt si ottengono le equazioni di
modulazione della soluzione oscillante. Tali equazioni, dette equazioni di Whitham, descrivono completamente l’evoluzione dei parametri lentamente variabili
della soluzione di propagazione.
Nel lavoro originale [1], Whitham mostrò l’applicazione della teoria per la
Korteweg-de Vries (KdV). Nella KdV l’evoluzione dei parametri della soluzione di propagazione è legata ad un sistema di tre PDE non lineari accoppiate.
Whitham dimostrò che tali equazioni possono essere diagonalizzate, tramite manipolazioni algebriche non banali, ottenendo un sistema di tre PDE iperboliche
quasi-lineari. Si sottolinea che in generale la diagonalizzazione dei sistemi di
PDE è limitata ai sistemi iperbolici 2 x 2 [2], quindi una classe di equazioni mol2
3
to limitata. Resta un problema aperto la soluzione dei sistemi con un numero
di equazioni maggiore di 2.
La risoluzione completa di un problema di DSW passa per la costruzione di
Gurevich e Pitaevskii [3] dove si impongono le condizioni iniziali alle equazioni di modulazione. Essi risolsero il problema di Riemann per la KdV, ovvero
la soluzione generale a partire da condizioni iniziali composte di due costanti
separate da una discontinuità finita. L’obiettivo della costruzione di GurevichPitaevskii è determinare la forma analitica delle velocità con cui si propaga lo
shock all’interno del mezzo. Tradizionalmente questo tipo di condizioni iniziali
è legata a problemi di gas dinamica e shock viscosi. L’applicazione del problema
di Riemann alla KdV unita con la costruzione di Gurevich-Pitaevskii conduce
alla completa risoluzione dello shock dispersivo.
Successivamente, la teoria di Whitham è stata estesa a soluzioni multifase
tramite l’uso della trasformata spettrale [4] e l’integrazione per gap finiti [5]. I
benefici di questo approccio includono la diagonalizzazione diretta delle equazioni di modulazione. Questa formulazione ha chiarito che la teoria di Whitham
è un metodo applicabile esclusivamente alla classe di equazioni che sono integrabile tramite il metodo della trasformata spettrale. Tale proprietà ha permesso
l’applicazione ad altre equazioni integrabili come la nonlinear Schrödinger (NLS)
[6].
Il limite della teoria di Whitham risiede nella sua scarsa applicabilità. Il
suo utilizzo è circoscritto ai modelli integrabili tramite trasformata spettrale,
mentre, a parte rari casi, i modelli fisici non possiedono questa proprietà. Tuttavia, nel 2005 [7], El ha introdotto un metodo, basato sulla teoria di Whitham
e la costruzione di Gurevich-Pitaevskii, per il calcolo delle velocità dei DSW
nei modelli non integrabili. Questo metodo utilizza alcune proprietà legate al
raccordo tra le equazioni di modulazione non diagonalizzate e il limite non dispersivo del modello, portando ad una nuova formulazione della costruzione di
Gurevich-Pitaevskii. Si vuole sottolineare che a differenza della teoria originale,
questo metodo non caratterizza completamente il profilo di modulazione, ma
risolve lo shock dispersivo senza conoscere la forma diagonale delle equazioni di
modulazione.
Nei capitoli che seguono verrà esposta in dettaglio la soluzione generale dello
shock dispersivo relativo alla KdV, in linea con gli articoli originari di Whitham
e Gurevich-Pitaevskii. Sulla base di questi risultati si introdurrà l’estensione
della teoria di Whitham ai modelli non integrabili, a partire dal lavoro di El
[7]. Si mostrerà che questa formulazione applicata alla KdV porta agli stessi
risultati ottenibili tramite la teoria di Whitham originale.
Una volta introdotta la teoria, si passera alla parte originale del lavoro di tesi.
In particolare, sarà presentata l’applicazione concreta della teoria di Whitham
per la nonlinear Schrödinger quintica (CQNLS), con una nonlinearità defocusing
per il termine cubico e una nonlinearità focusing nel termine quintico in regime
semiclassico (non linearità maggiore della dispersione). In aggiunta alla teoria di
Whitham, verrà esposto uno studio approfondito sul meccanismo di bistabilità
e instabilità che lega le soluzioni solitoniche della CQNLS.
Capitolo 2
Teoria di Whitham
La teoria di Whitham, nella sua formulazione originale, permette lo studio
analitico degli shock dispersivi, limitato nei modelli integrabili. Verrà presentata
l’applicazione della teoria per la Korteweg - de Vries (KdV), la più semplice
equazione che presenta la formazione di DSW.
2.1
Formazione di shock
La più semplice equazione di propagazione non lineare e non dispersiva che
presenta formazione di shock è l’equazione di Hopf
ut + uux = 0.
(2.1)
Una possibile risoluzione dell’equazione di Hopf avviene tramite il metodo
delle caratteristiche [2]. Prendendo una parametrizzazione del campo di velocità
u(x, t) tramite un parametro s
u = u(x(s), t(s)),
con un curva C
!
x = x(s)
t = t(s)
(2.2)
(2.3)
la derivata direzionale calcolata lungo la curva C é
du
du dt
du dx
=
+
.
ds
dt ds dx ds
(2.4)
Richiedendo che du/ds = 0 si ricava dal confronto con la (2.1)
dt
= 1,
ds
dx
= u.
ds
(2.5)
Si passa da un’equazione alle derivate parziali ad due equazioni ordinarie. Dalla (2.5), si ricava facilmente che t=s e la curva C si riduce a x=x(t).
L’equazione di Hopf espressa come derivata direzionale é
d
[u(x(t), t] = 0.
dt
4
(2.6)
5
2.1. FORMAZIONE DI SHOCK
La u [x(t), t] si mantiene costante lungo la curva caratteristica. Passando al
problema di Cauchy, prendendo come curva iniziale l’asse reale e il valore iniziale
della coordinata nel punto x=x∗
!
x(0) = x∗
(2.7)
u(x, 0) ≡ u0 = F (x)
Risolvendo il sistema di equazioni ordinarie (2.5) si ottiene
!
x(t) = x∗ + F (x)t
u(x) = F (x∗ )
(2.8)
data l’arbitrarietà di x∗ , si ha una famiglia di soluzione al variare della coordinata iniziale. Fissati x e t si può risalire alla u:
u(x, t) = F [x − u(x, t)t]
(2.9)
u!x, t"
1.0
0.8
0.6
0.4
0.2
x
!2
!1
1
2
3
Figura 2.1: Soluzione dell’equazione di Hopf (2.9) a diversi tempi, con F (x) =
1/(1 + x2 )
Si può riscrivere la soluzione nella forma:
x − ut = W (u),
(2.10)
dove W (u) è la funzione inversa della distribuzione iniziale F (x). Le curve
caratteristiche (2.10) hanno la forma di fasci di rette con coefficiente angolare
u(x, t).
Dopo un tempo tc la rette si intersecano e la soluzione diventa a piú valori
sviluppando una discontinuità. La dinamica descritta dalla (2.1) é caratterizzata dalla formazione di singolarità per t = tb , x = xb , in cui la soluzione formale
6
2.2. TEORIA DI WHITHAM: FORMULAZIONE GENERALE
2.0
1.5
1.0
0.5
5
10
15
20
Figura 2.2: Grafico delle rette caratteristiche dell’eq. di Hopf con F (x) =
1/(1 + x2 ). Notare l’intersezione che porta alla formazione dello shock.
diventa a tre valori e perde di significato fisico.
Il tempo in cui si forma la singolarità si ottiene derivando la (2.9)
du
= F " (1 − tux ) = 0
dx
→
tb = −
1
.
F"
(2.11)
Per studiare la dinamica oltre al singolarità, si includono nell’equazione effetti
di ordine superiore, tali da trasformare l’equazione di Hopf in:
• Burgers: ut + uux + νuxx = 0
correzione legata alla viscosità
0<ν$1
• Korteweg - de Vries: ut + uux + #uxxx = 0
correzione legata alla dispersione
0<#$1
La risoluzione delle due equazioni conduce a soluzioni e situazioni fisiche profondamente diverse. Mentre nei modelli viscosi si osserva uno smorzamento dello
shock, i modelli dispersivi presentano la formazione di una regione oscillatoria
che generalmente termina in un solitone.
2.2
Teoria di Whitham: formulazione generale
Un sistema non lineare che presenta una soluzione di propagazione si possono distinguere due scale di evoluzione: la scala di oscillazione veloce (θ = x − ct)
e la scala lenta di modulazione. Si cerca una soluzione per onda cnoidale dipendente da alcuni parametri lentamente variabili. Le equazioni che descrivono
l’evoluzione di questi parametri sono dette equazioni di modulazione (o di Whitham) e si trovano mediando le leggi di conservazione sulla scala di oscillazione
veloce [8].
7
2.2. TEORIA DI WHITHAM: FORMULAZIONE GENERALE
Si consideri una equazione generica
Φ(u, ut , ux , uxx , uxt , ..) = 0,
(2.12)
supponendo che la variabile u(x, t) dipenda da x e t solo nella combinazione
u = u(θ),
con θ = x − ct,
(2.13)
l’ Eq. (2.12) si riduce nella forma
u2θ = F (u, c, Ai ),
(2.14)
dove c rappresenta la velocità di propagazione e le Ai di integrazione introdotte
nel processo di integrazione della (2.12) (es. energia, quantità di moto). Le
soluzioni periodiche corrispondono al moto di u tra due zeri della F (u, c, Ai ),
u1 (c, Ai ) e u2 (c, Ai ), con u1 < u2 , tali che, in accoro con la (2.14), la funzione
F sia positiva nell’intervallo u1 < u < u2 . Siano θ1 e θ2 i valori di θ per cui u
assume valore u1 e u2
u(θ1 ; c, Ai ) = u1 ,
u(θ2 ; c, Ai ) = u2 ,
(2.15)
integrando per quadratura la (2.14) si ricava la lunghezza d’onda L
L=2
"
θ2
dθ = 2
θ1
"
u2
u1
du
=2
uθ
"
u1
Corrispondentemente si ottengono le relazioni:
k = k(c, Ai ) =
2π
,
L(c, Ai )
u2
du
#
.
F (u, c, Ai )
ω = ω(c, Ai ) = ck(c, Ai ),
(2.16)
(2.17)
dove i parametri c e Ai sono funzioni lente di x e t e hanno una variazione lenta
in un periodo di 2π/ω. Si definisce l’operazione di media come
f=
1
2∆
"
x+∆
f (x" , t)dx" .
(2.18)
x−∆
supponendo che l’intervallo 2∆ sia maggiore della lunghezza d’onda L.
L’ ipotesi di Whitham è che la mediazione delle leggi di conservazione permetta di ricavare le equazioni di modulazione dei parametri lentamente variabili
della soluzione (2.14). Quindi, nota una legge di conservazione
ρt + Jx = 0,
(2.19)
si applica la mediazione definita nella (2.18). Per la ρ si ha
e per la J
$ %
J =
x
" x+∆
$ %
$ %
1
∂
ρ =
ρ(x" , t)dx" = ρ ,
2∆ x−∆ ∂t
t
t
1
2∆
%
J(x" + ∆, t) − J(x" − ∆, t) =
$ %
& x+∆
∂ 1
"
"
= ∂x
.
2∆ x−∆ J(x , t)dx = J
& x+∆
∂
(x" , t)dx"
x−∆ ∂x
=
1
2∆
$
x
(2.20)
(2.21)
8
2.3. TEORIA DI WHITHAM NELLA KDV
Le leggi di conservazione mediate
ρt + J x = 0.
(2.22)
danno la variazione dei parametri sulla scala lenta di modulazione. Si devono
avere un numero di leggi di conservazione pari al numero di parametri che descrivono il sistema. Si ha in questo modo un sistema, generalmente complicato,
la cui risoluzione porta ad una descrizione completa dell’onda cnoidale.
2.3
Teoria di Whitham nella KdV
Si vuole applicare la teoria di Whitham alla più semplice equazione dispersiva
che presenta lo shock, la KdV. Quindi seguendo il procedimento sopra descritto:
trovare una soluzione del tipo (2.14), definire l’operazione di media sulla scala
veloce e ricavare le equazioni di modulazione dalle leggi di conservazione.
2.3.1
Onde Cnoidali
Il primo passo per l’applicazione della teoria di Whitham è il calcolo delle
onde cnoidali della KdV. Si consideri la KdV unidimensionale, per tempi in cui
il termine dispersivo è confrontabile con il termine non lineare (# ≈ 1)
ut + uxxx + uux = 0.
(2.23)
Sostituendo la soluzione di propagazione u(θ), con θ = x − ct, si ha
−cuθ + uθθθ + uuθ = 0
(2.24)
e procedendo con una integrazione rispetto a θ
−cu + uθθ +
u2
= B,
2
(2.25)
con B costante di integrazione. Si riconosce nella (2.25) un sistema Hamiltoniano. Definendo uθ = p, si ha che le variabili canoniche (p, u) devono soddisfare
le equazioni di Hamilton

2
 pθ = −∂H/∂u = B + cu − u2
(2.26)

uθ = ∂H/∂p = p
da cui si ottiene l’Hamiltoniana H:
H=
p2
cu2
u3
− Bu −
+ .
2
2
6
(2.27)
Per convenienza si definisce una nuova costante A ≡ −H e l’equazione (2.27) si
riscrive
u2θ = −G(u).
(2.28)
dove G(u) è un polinomio cubico
G(u) ≡ 2A − 2Bu − cu2 +
u3
.
3
(2.29)
9
2.3. TEORIA DI WHITHAM NELLA KDV
Indichiamo con u1 ,u2 e u3 i tre zeri reali di G(u) con u3 > u2 > u1 . Con questa
posizione la (4.53) diventa
G(u) =
1
(u − u1 )(u − u2 )(u − u3 ),
3
(2.30)
Confrontando la (4.53) con la (4.143) si ricavano le espressioni delle costanti del
moto in funzione dei tre zeri di G(u)
c=
1
(u1 + u2 + u3 )
3
(2.31)
1
A = − (u1 u2 u3 )
6
1
B = − (u1 u2 + u2 u3 + u1 u3 ).
6
La soluzione della (2.28) si ottiene per quadratura;
du
dθ = #
.
−G(u)
(2.32)
(2.33)
(2.34)
La (2.34) ammette soluzioni reali solo se il radicando è maggiore di zero. Uno
studio della G(u) rivela che la richiesta è soddisfatta per valori di u tali che
u2 < u < u 3 .
0.2
0.15
0.1
G(u)
0.05
0
−0.05
−0.1
−0.15
−0.2
−0.5
0
0.5
u
Figura 2.3: Andamento del polinomio G(u) per tre zeri generici.
1
10
2.3. TEORIA DI WHITHAM NELLA KDV
La (2.34) si integra tra u2 = u(θ2 ) e u3 = u(θ3 ) ottenendo
" u3
du
#
(θ3 − θ2 ) =
.
−G(u)
u2
(2.35)
Le condizioni periodicità dell’onda cnoidale impone:
u(θ + L) = u(θ),
dove si definisce L, il periodo, come
"
L≡2
u3
u2
#
(2.36)
du
−G(u)
.
Il vettore d’onda k corrispondente è
* " u3
+−1
du
2π
#
k=
=π
,
L
−G(u)
u2
(2.37)
(2.38)
e dalla (2.17) si definisce la relazione
ω=
1
(u1 + u2 + u3 )k.
3
(2.39)
Data la fattorizzazione del polinomio G(u) nella (4.143), la soluzione della
(2.35) è esprimibile tramite le funzioni ellittiche [9]. Operando la sostituzione
u = u3 − (u3 − u2 ) sin2 φ"
(2.40)
il polinomio G(u) diventa:
%$
$u − u %
%
$
%2 1 ,$
3
2
, (2.41)
G(u) =
u3 − u1 1 −
sin2 φ" sin2 φ cos2 φ" u3 − u2
3
u3 − u1
con il differenziale:
du = 2(u2 − u3 ) cos φ" sin φ" dφ"
Sostituendo nella (2.35) si ottiene:
√
" φ
2 3
dφ"
#
θ=√
u3 − u1 0
1 − m sin φ" 2
(2.42)
(2.43)
dove m, modulo, assume valore 0 ≤ m ≤ 1 ed è definito
m=
u3 − u2
.
u3 − u1
(2.44)
La soluzione espressa in funzione dell’integrale ellittico di Jacobi K(φ, m) ha la
forma
.
u3 − u1
K(φ, m) ≡
θ.
(2.45)
3
Il vettore d’onda k diventa
1
π (u3 − u1 ) 2
k= √
2 3 K(φ, m)
(2.46)
11
2.3. TEORIA DI WHITHAM NELLA KDV
Invertendo la (2.40) si ricava l’espressione della soluzione ellittica per la KdV
/.
0
u3 − u1
2
u(x, t) = u3 − asn
(x − ct); m ,
(2.47)
3
dove a è l’ampiezza e pari a
a = u3 − u2 .
(2.48)
sn2 (φ; m) = sin2 φ,
(2.49)
La soluzione ottenuta dipende dal parametro m, è utile studiare il comportamento dell’onda cnoidale nei due limiti imposti dal modulo.
Se m → 0, si ha che u3 ≈ u2 e l’onda cnoidale (2.47) si trasforma in una funzione
trigonometrica. Il seno ellittico di Jacobi diventa
L’integrale ellittico del primo tipo K(φ, m), assume valore K(φ, 0) = π2 e il
vettore d’onda k diventa:
.
u3 − u1
k0 ≡ k(m = 0) =
.
(2.50)
3
Sostituendo la (2.50) nella (2.31), si ricava l’espressione della velocità c in questo
limite
c0 = u3 − k02 .
(2.51)
Quindi, l’espressione dell’onda cnoidale nel limite m = 0 è
,
u(x, t) = u3 − a sin2 k0 (x − c0 t) ,
(2.52)
Dalla (2.39) si definisce la relazione di dispersione a cui è soggetta la (2.52)
ω0 = co k0 = u3 k0 − k03 .
(2.53)
La stessa relazione è ricavabile dalla linearizzazione della KdV per piccole
ampiezze intorno ad un valore di fondo u = u3 .
Si può concludere, quindi, che il limite m = 0 è il limite lineare dell’onda
cnoidale.
Nel limite m → 1, u1 ≈ u2 e l’onda cnoidale si trasforma in una funzione
iperbolica
1
cn(θ, 1) =
,
(2.54)
cosh θ
L’ampiezza a diventa
as = u3 − u1 ,
(2.55)
il numero d’onda k tende a 0 e il periodo L tende a infinito. Dalle relazioni sulle
funzioni di Jacobi, usando l’identità
sn(θ, m)2 + cn(θ, m)2 = 1,
(2.56)
e la (2.55), sostituendo nella (2.47), si ottiene il solitone della KdV
u(x, t) = u1 +
a
# s
,
cosh ( as /3(x − cs t))
2
(2.57)
12
2.3. TEORIA DI WHITHAM NELLA KDV
dove, la velocità di traslazione del solitone, espressa in funzione di as ,
cs = u1 +
as
.
3
(2.58)
Quindi il limite m = 1 rappresenta il limite solitonico dell’onda cnoidale.
Nelle equazioni trovate (2.57) e (2.52) i termini costanti u1 e u3 generalmente si
eliminano ponendo la condizione naturale u(x, t) = 0 nei limiti x → ±∞. Nella
teoria di modulazione i termini costanti non possono essere eliminati perché i tre
zeri u1 , u2 e u3 sono funzioni lente di x e t e contribuiscono nella regolarizzazione
dello shock.
2.3.2
Equazioni di modulazione
Il secondo passo della teoria è determinare le equazioni di evoluzione delle
costanti del moto A, B in aggiunta alla velocità c. Le equazioni di modulazione
si ricavano a partire dalle leggi di conservazione. Il numero di leggi necessario è
pari al numero di parametri che definiscono l’onda cnoidale e la mediazione su
un periodo (2.37) descrive la dinamica sulla scala lenta dell’onda cnoidale. La
riduzione della KdV necessita due costanti, A e B, che aggiunte alla velocità
c diventano i tre parametri dell’onda cnoidale. Tali parametri sono legati ai
tre zeri del polinomio G(u) dalla (2.31)-(2.33), possiamo quindi formulare il
problema in funzione delle variabili u1 , u2 e u3 .
La prima si ottiene semplicemente dalla ristruttura della KdV ed esprime la
conservazione della massa:
%
$ %
$ u2
u +
+ uxx = 0.
(2.59)
2
t
x
Una seconda legge, legata alla conservazione della quantità u2 /2, esprime la
conservazione dell’energia,
$ u2 %
2
t
+
$ u3
6
+ uuxx −
u2x %
= 0.
2 x
(2.60)
La terza legge di conservazione, legata all’invarianza dell’ Hamiltoniana, è ottenuta moltiplicando la (2.23) per u2 .
$ u3
6
−
$ u4
%
u2x %
u2
u2
+
+ uxx + xx + ut ux = 0.
2 t
8
2
2
x
(2.61)
Nella KdV è possibile esprimere tutte le leggi di conservazione in funzione
di poche grandezze mediate, semplificando notevolmente le tre leggi di conservazione. Tale proprietà deriva direttamente dalla integrabilità della KdV.
In virtù della soluzione u(x − ct), si ha che ∂θ = ∂x . Quindi, dalla (2.59)
sostituendo al posto di uθθ la (2.25) si ottiene per la prima legge:
$ %
$
%
u + B + cu = 0,
(2.62)
t
x
Per la seconda legge (2.60) sostituendo prima la (2.25) e successivamente la
(2.28) si ricava:
$ u2 %
$ cu2
%
+
+A .
(2.63)
2 t
2
x
13
2.3. TEORIA DI WHITHAM NELLA KDV
Per la terza sostituendo al primo termine la (2.27) e la (2.23) al posto del termine
ut ux si ottiene:
$
$ B2
%
%
cu2
c2 u2
− u2θ +
+ Bcu +
− cu2θ = 0
2
2
2
t
x
− A + Bu +
(2.64)
Definiamo la media di una osservabile F (u) su un periodo L come:
1
F (u) =
L
"
L
0
k
dθF (θ; u1 , u2 , u3 ) =
2π
"
2π
k
0
F (u)
du #
.
−G(u)
(2.65)
Il calcolo di un integrale, del tipo definito nella (2.65), è generalmente complicato. É riportato un esempio dell’applicazione per il calcolo di u. Per definizione
si ha:
"
k u3
u
u=
du #
(2.66)
π u2
−G(u)
operando la sostituzione:
u = u3 − (u3 − u2 ) sin2 φ = u1 + (u3 − u1 )(1 − m sin2 φ)
(2.67)
otteniamo
u=
k
π
"
u3
u2
,
1
du #
u1 + (u3 − u1 )(1 − m sin2 φ) .
−G(u)
(2.68)
Per il primo termine u1 è costante e l’integrale restante è pari a uno per definizione. Per il secondo termine invece usando le espressioni (2.41) e (2.42)
otteniamo:
√
"
1
4 3k (u3 − u1 ) φ
√
u = u1 +
dφ(1 − m sin2 φ) 2
(2.69)
π
u3 − u1 0
Sostituendo la (2.46) otteniamo l’espressione:
u = u1 + 2(u3 − u1 )
E(m)
.
K(m)
(2.70)
& #
dove E(m) = dφ 1 − m sin2 φ, integrale ellittico del secondo tipo.
Applicando l’operazione di media definite nella (2.65) e in virtù della linearità
di tale operazione il sistema mediato diventa
$
%
 $ %

u + B + cu = 0


t
x








 $
%
$
%
1 2
1
2
u
+
A
+
cu
=0
2
2

t
x








$
%
$


 − A + Bu + c u2 − u2 + B 2 + Bcu +
θ
2
2
t
(2.71)
c2 2
2 u
− cu2θ
%
x
=0
14
2.3. TEORIA DI WHITHAM NELLA KDV
dove le quantità mediate da calcolare si sono ridotte a: 1, u, u2 e u2θ .
Seguendo il procedimento originario proposto da Whitham [8] si esprimono le
quantità mediate in funzione di derivate rispetto alle costanti della funzione
ellittica W = W (A, B, c) definita come
1
W (A, B, c) = −
2π
"
0
L
#
−G(u)du.
(2.72)
La G(u) è definita nella (4.53), A e B sono le costanti note dalla riduzione e c
è la velocità di propagazione. Le quantità medie sono
" L
k
1
1=
du #
= kWA ,
2π 0
−G(u)
" L
k
u
u=
du #
= −kWB ,
2π 0
−G(u)
" L
1 2
k
u2
u =
du #
= −kWc ,
2
2π 0
2 −G(u)
" L
k
u2
u2θ =
du # θ
= −kW.
2π 0
−G(u)
(2.73)
(2.74)
(2.75)
(2.76)
Per il momento si mantengono impliciti gli integrali (2.73), (2.76). Sostituendo nelle equazioni (2.71), si ottiene il sistema espresso in termini di
W (A, B, c).
%
$
$
%
kWB + ckWB − B = 0,
(2.77)
t
x
$
%
$
%
kWc + ckWc − A = 0,
(2.78)
t
,
x
- ,
1 k(AWA +BWB +cWC −W ) + A(kWA )t +B(kWB )+c(kWc )−W kt − B 2 = 0
2
t
x
(2.79)
La terza equazione può essere semplificata usando alcune proprietà di W .
Analizzando il primo membro
,
k(AWA + BWB + cWC − W ) ,
(2.80)
t
le derivate di W rispetto a x e t sono
∂t W (A, B, c) = WA At + WB Bx + Wc ct
∂x W (A, B, c) = WA Ax + WB Bx + Wc cx .
(2.81)
Derivando le costanti A, B e c si ottiene l’identità (2.81), per i restanti termini
abbiamo:
,
A(kWA )t + B(kWB ) + c(kWc ) − W kt
(2.82)
Per il secondo membro invece, aggiungendo e sottraendo (cA)x e poi sostituendo le funzioni in termini dell’integrale W (A, B, c) si ha:
,
1
kc(AWA + BWB + cWC − W ) − B 2 − Ac = 0.
2
x
(2.83)
15
2.3. TEORIA DI WHITHAM NELLA KDV
operando allo stesso modo si ricava:
,
A(ckWA )x + B(ckWB )x + c(ckWc )x − W (ck)x − BBx − (cA)x .
Quindi, sommando la (2.82) e la (2.84)
,
A (kWA )t + (ckWA − c)x +
,
+B (kWB )t + (ckWB − B)x +
,
+c (kWc )t + (ckWc − A)x −
,
−W (k)t + (ck)x = 0
(2.84)
(2.85)
Il primo termine è nullo per la (2.73), il secondo e il terzo invece si annullano
per la (2.77) e la (2.78). Si ottiene, in definitiva, una legge di conservazione
generale
(k)t + (ω)x = 0.
(2.86)
con ω = ck, che esprime la conservazione del numero d’onda. Svolgendo le
derivate:
+
*
∂
∂
+c
k = −kcx ,
(2.87)
kt + cx k + ckx = 0 →
∂t
∂x
Inoltre, dato che vale la (2.73), si ricavano le equazioni di modulazione:
Dt WA = WA cx ,
Dt WB = WA Bx ,
Dt Wc = WA Ax .
(2.88)
dove Dt = (∂/∂t + c∂/∂x).
Si vuole esprimere il sistema (2.88) in funzione degli zeri della G(u). Dalle
relazioni (2.31), (2.33) e (2.32) si ricava:
cx =
1,
u1 + u2 + u3 ,
3
x
(2.89)
1,
(u1 )x (u2 + u3 ) + (u2 )x (u3 + u1 ) + (u3 )x (u1 + u2 ) ,
(2.90)
6
1,
Ax = − (u1 )x u2 u3 + (u2 )x u2 u3 + (u3 )x u2 u1 .
(2.91)
6
Dato che A, B e c sono funzioni u1 , u2 e u3 si ha che W = W (u1 , u2 , u3 ) e
il sistema (2.88) diventa
Bx = −
WA,u1 Dt u1 + WA,u2 Dt u2 + WA,u3 Dt u3 =
$
%
= 13 WA u1 + u2 + u3 ,
(2.92)
x
=
− 61 WA
W
D u + WB,u2 Dt u2 + WB,u3 Dt u3 =
$ B,u1 t 1
%
(u1 )x (u2 + u3 ) + (u2 )x (u3 + u1 ) + (u3 )x (u1 + u2 ) ,
x
(2.93)
16
2.3. TEORIA DI WHITHAM NELLA KDV
Wc,u1 Dt u1 + Wc,u2 Dt u2 + Wc,u3 Dt u3 =
$
%
= − 16 WA (u1 )x u2 u3 + (u2 )x u2 u3 + (u3 )x u2 u1 .
(2.94)
x
dove:
$
$
WA
WB
$
Wc
%
%
%
1
4π
1
4π
u1
=−
u1
1
=−
4π
e analogamente per u2 e u3 .
2.3.3
u1
=
"
"
"
du
#
(u − u1 ) −G(u)
(2.95)
du
#
.
(u − u1 ) −G(u)
(2.97)
du
#
,
(u − u1 ) −G(u)
(2.96)
Diagonalizzazione delle leggi di conservazione
La diagonalizzazione del sistema si esegue tramite la ricerca degli invarianti
di Riemann. Si consideri la somma delle (2.92)-(2.94) moltiplicate rispettivamente per tre parametri arbitrari p, q e r
$
%
(2.98a)
p(WA )u1 + q(WB )u1 + r(Wc )u1 Dt u1 +
$
%
+ p(WA )u2 + q(WB )u2 + r(Wc )u2 Dt u2 +
(2.98b)
$
%
+ p(WA )u3 + q(WB )u3 + r(Wc )u3 Dt u3 =
(2.98c)
%
2
$
1
= WA 2p u1 + u2 + u3 −
(2.98d)
6
x
−q [(u1 )x (u2 + u3 ) − (u2 )x (u3 + u1 ) − q(u3 )x (u1 + u2 )]
(2.98e)
3
−r [(u1 )x u2 u3 − r(u2 )x u2 u3 − r(u3 )x u2 u1 ]
(2.98f)
Focalizzando l’attenzione sulle eq (2.98(d-f)), il termine a destra dell’uguale si
scrive
,
1
(u1 )x (2p − q(u2 + u3 ) − ru2 u3 ) + (u1 )x (2p − q(u1 + u3 ) − ru1 u3 )+
W
A
6
+(u1 )x (2p − q(u2 + u1 ) − ru2 u1 ).
(2.99)
Scegliendo i parametri p, q e r tali che u3,x = 0 e u2,x = u1,x si ottiene
!
2p − q(u2 + u1 ) − ru2 u1 = 0
(2.100)
2p − q(u1 + u3 ) − ru1 u3 = 2p − q(u2 + u3 ) − ru2 u3
il cui risultato da:

 p = u1 u2 − u3 u2 − u3 u2
q = 2u3

r = −2
(2.101)
dove r è stato scelto in modo arbitrario. In questo modo le (2.98(d-f)) diventano
1
∂
− (u3 − u2 )(u3 − u1 )WA (u1 + u2 ).
3
∂x
(2.102)
2.3. TEORIA DI WHITHAM NELLA KDV
17
Per le (2.98(a-c)) si ha per ogni derivata totale Dt un somma dei tre funzionali WA , WB e Wc . Si consideri il termine a moltiplicare Dt U3 . Sostituendo
le espressioni in funzione di W (A, B, c) e dei parametri p, q e r definiti nella
(2.101) si ottiene:
1
4π
"
du
#
(u3 u2 + u3 u1 − u2 u2 − 2uu3 + u2 ) =
(u − u3 ) −G(u)
4
"
1
d ,
(u − u2 )(u − u1 ) =
2
du = 0
4π
du
−(u − u3 )
Per Dt u2 si ha
"
1
du
#
(u3 u2 + u3 u2 − u1 u2 − 2uu3 + u2 ) ≡ k1 ,
4π
(u − u2 ) −G(u)
mentre per Dt u3
"
1
du
#
(u3 u2 + u3 u2 − u1 u2 − 2uu3 + u2 ) ≡ k2 ,
4π
(u − u1 ) −G(u)
(2.103)
(2.104)
(2.105)
(2.106)
dalle espressioni precedenti si verifica che
k1 − k2 = 0.
(2.107)
É possibile riscrivere k1 e k2 in modo tale da estrarre un differenziale esatto
simile alla (2.104)
1
+ 4π
e
1
+ 4π
&
&
−u2 )
k1 = − (u32π
du
√
(u2
(u−u2 ) −G(u)
u−u1
u−u2
√
du
+
−G(u)
(2.108)
du
+
−G(u)
(2.109)
− 2uu2 + u3 u2 + u1 u2 − u3 u1 )
−u1 )
k2 = − (u32π
du
√
(u2
(u−u1 ) −G(u)
&
&
u−u2
u−u1
√
− 2uu1 + u3 u1 + u1 u2 − u3 u2 ),
ovviamente il secondo termine nelle due espressioni nullo perché si effettua una
integrazione su un percorso chiuso. I termini restanti, aggiungendo e sottraendo
la (2.73), possono essere riscritti come
k1 = −(u3 − u2 )WA − 2(u3 − u2 )(u2 − u1 )WA,u2
k2 = −(u3 − u1 )WA − 2(u3 − u1 )(u2 − u1 )WA,u1
Dato che k1 − k2 = 0 si ha:
$
%
WA = −2 (u3 − u2 )WA,u2 + (u3 − u1 )WA,u1 ,
(2.110)
(2.111)
quindi, sostituendo in k1 o in k2
k1 = k2 = −2(u3 − u2 )(u3 − u1 )(WA,u1 + WA,u2 ),
(2.112)
18
2.3. TEORIA DI WHITHAM NELLA KDV
Dalla definizione di G(u) e WA si nota che vale l’identità
WA,u1 + WA,u2 + WA,u3 ≡ 0
(2.113)
e la (2.112) si riscrive come:
2(u3 − u2 )(u3 − u1 )WA,u3
(2.114)
Eguagliando la (2.102) e la (2.114) otteniamo l’espressione ridotta della
(2.88)
1
∂
2(u3 − u2 )(u3 − u1 )WA,u3 Dt (u1 + u2 ) = − (u3 − u2 )(u3 − u1 )WA (u1 + u2 ),
3
∂x
(2.115)
quindi,
1 WA ∂
Dt (u1 + u2 ) +
(u1 + u2 ) = 0
(2.116)
6 WA,u3 ∂x
Allo stesso modo permutazioni cicliche di u1 , u2 e u3 danno le equazioni:
Dt (u1 + u3 ) +
1 WA ∂
(u1 + u3 ),
6 WA,u2 ∂x
(2.117)
Dt (u3 + u2 ) +
1 WA ∂
(u3 + u2 )
6 WA,u1 ∂x
(2.118)
Definiamo ora i seguenti parametri, noti come invarianti di Riemann
u1 + u2
,
2
u1 + u3
β2 =
,
2
u2 + u3
β3 =
.
2
Invertendo le relazioni (2.116)-(2.117) in funzione degli invarianti βi
β1 =
(2.119)
(2.120)
(2.121)
si ricava
u1 = β1 + β2 − β3 ,
(2.122)
u2 = β1 − β2 + β3 ,
(2.123)
u3 = −β1 + β2 + β3 .
(2.124)
WA,β1 = WA,u1 + WA,u2 − WA,u3
(2.125)
WA,β1 = −2WA,u3 .
(2.126)
Quindi le derivate di WA = WA (u1 , u2 , u3 ) espresse in funzione degli invarianti diventano
e usando la (2.113), si ottiene
Allo stesso modo si ricava per le altre due espressioni
WA,β3 = −2WA,u1 ,
WA,β2 = −2WA,u2 .
Dato che, WA = 1/k = L/2π, si verifica che
(2.127)
19
2.3. TEORIA DI WHITHAM NELLA KDV
WA
2WA
2
2k
=
=−
=
.
WA,u1
WA,β3
∂(ln WA )/∂β3
kβ3
(2.128)
Allo stesso modo si ottengono le relazioni per WA,u1 eWA,u2 . Sostituendo le
relazioni trovate nelle Eq (2.116)e (2.117) e espandendo Dt si ricava
∂
∂
βi + Vi (β1 , β2 , β3 ) βi
∂t
∂x
con
Vi =
i = 1, 2, 3.
1
1
(β1 + β2 + β3 ) −
,
3
3∂βi (ln L)
dove m e L si ricavano dalla (2.46), e pari a
√
2 6
L= √
K(φ, m),
β3 − β1
m=
β2 − β1
.
β3 − β1
(2.129)
(2.130)
(2.131)
Con K(φ, m) integrale ellittico del primo tipo. Per il calcolo del secondo termine
della Vi si ha
∂K(φ, m)
∂K(φ, m) ∂m
=
(2.132)
∂βi
∂m
∂βi
e la derivate di K(φ, m) rispetto a m è
dK(m)
E(m) − (1 − m)(K(m))
=
;
dm
2m(1 − m)
(2.133)
con E(m) integrale ellittico del secondo tipo. Applicando le condizioni (2.132)(2.133) si ottengono le velocità Vi
V1 =
1
2 (β3 − β1 )(1 − m)K(m)
(β1 + β2 + β3 ) +
3
3
E(m) − K(m)
(2.134)
V2 =
1
2 (β2 − β1 )(1 − m)K(m)
(β1 + β2 + β3 ) −
3
3 E(m) − (1 − m)K(m)
(2.135)
V3 =
2 (β3 − β1 )(1 − m)K(m)
1
(β1 + β2 + β3 ) +
3
3
E(m)
(2.136)
dette velocità caratteristiche. Utilizzando le (2.122)-(2.124) è possibile esprimere
le onde cnoidali per la KdV in funzione degli invarianti di Riemann
$
%2 β + β − β
,#
-3
1
2
3
u = β3 − β1
− 2m sn2
β3 − β1 (x − ct), m
β3 − β1
e conseguentemente la velocità c, il modulo m e il vettore d’onda k
√
π β3 − β1
β2 − β1
k=√
,
m=
,
β3 − β1
6 K(φ, m)
c=
1
(β1 + β2 + β3 ),
3
(2.137)
(2.138)
(2.139)
20
2.4. COSTRUZIONE DI GUREVICH-PITAEVSKII
2.4
Costruzione di Gurevich-Pitaevskii
La costruzione di Gurevich-Pitaevskii fu introdotta originariamente nel 1974,
ponendosi il problema di raccordare la condizioni iniziali della KdV con le equazioni di Whitham. Questa costruzione ha carattere più generale ed è possibile
l’estensione a problemi dinamici più complicati (NLS).
Prendendo piccoli effetti di dispersione hanno mostrato che dopo un tempo finito, confrontabile con il tempo di rottura tb (2.11), si forma un fronte oscillatorio
seguito dalla formazione di un solitone [3].
Più precisamente, si consideri la KdV unidimensionale
ut + uux + uxxx = 0,
(2.140)
passando il sistema alla scala spaziale lenta X = #x con # $ 1 si ha che, durante
l’evoluzione iniziale, |ux | ∼ #, |uxxx | ∼ #3 . Quindi |uxxx | $ |ux | e nella prima
fase dell’evoluzione si può trascurare il termine dispersivo nella (2.140). In
questo limite, detto limite semiclassico, l’evoluzione è descritta dall’equazione
di Hopf con u(x, t) = f (x, t) e fxxx = 0,
f (x, t)t + f (x, t)f (x, t)x = 0.
(2.141)
L’evoluzione dell’equazione di Hopf porta alla formazione di singolarità e ad
una regione a più valori, come descritto nel paragrafo $2.1. L’ idea generale
della costruzione di Gurevich-Pitaevskii (GP) è di usare la teoria di Whitham
quando i termini dispersivi non sono più trascurabili e la soluzione acquista
carattere oscillatorio. In questa modellizzazione si può dividere il piano (x, t) in
tre domini, delimitati da due curve x− (t) e x+ (t):
(−∞, x− (t)) ∪ [x− (t), x+ (t)] ∪ (x+ (t), +∞)
t
Dispersionless limit
(2.142)
x (t )
Whitham
equations
x (t )
Dispersionless limit
x
Figura 2.4: Grafico qualitativo dei limiti di validità delle equazioni
2.4. COSTRUZIONE DI GUREVICH-PITAEVSKII
21
Al di fuori dell’intervallo compreso tra x− (T ) e x+ (T ) la soluzione è regolata
dall’equazione di Hopf, limite non dispersivo della KdV. La dinamica all’interno
dell’intervallo è localmente periodica e descritta dalle equazioni di Whitham.
Sui bordi, per x = x− (t), detto ”trailing edge”, l’onda di propagazione è lineare
con modulo m = 0 e ampiezza a = 0. Per x = x+ (t), detto ”leading edge”,
l’onda di propagazione diventa un solitone con m = 1 e numero d’onda k = 0.
Per verificare la validità di tale costruzione è necessario calcolare i limiti
armonico e solitonico della equazioni di modulazione. Se m = 0, dalla (2.138),
si ha β2 = β1 . Il calcolo diretto delle velocità nelle equazioni (2.134), (2.136) da
V3 = β3 ,
V2 = V1 = 2β1 − β3 .
(2.143)
Le equazioni di Whitham si riducono a
∂
∂
β3 + β3 β3 = 0,
∂t
∂x
∂
∂
β1 + (2β1 − β3 ) β1 = 0.
∂t
∂x
(2.144)
Se m = 1, si ha β3 = β2 . Le velocità diventano
V1 = β1 ,
V2 = V3 =
1
(β1 + 2β3 ),
3
(2.145)
e le equazioni di Whitham
∂
∂
β1 + β1 β1 = 0,
∂t
∂x
∂
1
∂
β3 + (β1 + 2β3 ) β3 = 0.
∂t
3
∂x
(2.146)
Una delle equazioni di Whitham diventa sempre l’equazione di Hopf (2.1), il limite non dispersivo della KdV. Usando la definizione degli invarianti di Riemann
si verifica facilmente, che le velocità caratteristiche multiple (2.143) e (2.144)
coincidono con la velocità di gruppo ricavabile dalla (2.53) nel limite armonico e
con cs (2.58) nel limite solitonico. Se una equazione diventa l’equazione di Hopf,
le restanti due collassano in una nuova equazione la cui velocità corrisponde o
alla velocità di gruppo del limite armonico o alla velocità del solitone KdV nel
limite iperbolico.
In conclusione, la costruzione di GP mostra una struttura pienamente compatibile con le equazioni di Whitham e tale compatibilità deve essere rispecchiata nelle condizioni di raccordo nei punti x = x± (t) tra la soluzione interna,
equazioni di Whitham, e la soluzione esterna, equazione di Hopf:
!
x = x− (T ) → β2 = β1 β3 = f (x, t),
(2.147)
x = x+ (T ) → β2 = β1 β1 = f (x, t),
dove f (x, t) è la soluzione dell’equazione di Hopf. Si sottolinea che i bordi
x± (t) non sono noti e la loro determinazione è parte stessa del problema.
Le velocità che descrivono la variazione dei bordi x− (t) e x+ (t) sono legate alle
velocità caratteristiche associate al collasso di due invarianti di Riemann nei due
limiti. In definitiva, i bordi sono descritti dall’equazione differenziale
dx−
= V2 (β1 , β1 , β3 )x=x− = V1 (β1 , β1 , β3 )x=x−
dt
dx+
= V2 (β1 , β3 , β3 )x=x+ = V3 (β1 , β3 , β3 )x=x+ .
dt
(2.148)
(2.149)
22
2.4. COSTRUZIONE DI GUREVICH-PITAEVSKII
Data l’invarianza del sistema per trasformazioni di scala x → Cx e t → Ct,
le equazioni di modulazione possono essere espresse in funzione della variabile
auto similare
x
ζ=
(2.150)
t
e ciò implica che
βi = βi (ζ)
(2.151)
Sostituendo la (2.151) nelle equazioni di Whitham (2.129) si ottiene
dβi
(Vi − ζ) = 0,
dζ
i = 1, 2, 3.
(2.152)
Questa costruzione è una costruzione generale. Il raccordo ai bordi tra la soluzione interna e la soluzione esterna è strettamente legato alle differenti condizioni
al contorno che è possibile imporre al sistema.
Prendendo come condizioni iniziali una discontinuità finita
u0 (−∞) = u− ,
u0 (+∞) = u+
con u− > u+ ,
(2.153)
Per queste condizioni, la soluzione auto similare della KdV nel limite non
dispersivo è:
 −
x < t,
 u
x/t,
0<x<t
u(x, t) =
(2.154)
 +
u ,
x > t,
che assume valore multiplo per x− (t) < x < x+ (t), come riportato in figura
(2.5). In accordo con la costruzione di Gurevich e Pitaevskii nella regione a va-
1
u=u−
0.8
u(x)
0.6
0.4
0.2
u=u+
0
x−(t)
x+(t)
−0.2
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
x
Figura 2.5: Grafico della soluzione della KdV nel limite non dispersivo per
u− = 1 e u+ = 0.
23
2.4. COSTRUZIONE DI GUREVICH-PITAEVSKII
lore multiplo si forma l’onda non lineare modulata, la cui evoluzione è governata
dalle equazioni di Whitham (vedi figura (2.6)). Gli invarianti di Riemann devono soddisfare il limite non dispersivo nella KdV, quindi, le condizioni (2.147)
diventano:
!
s = s− → β2 = β1 β3 = u− ,
(2.155)
s = s+ → β2 = β3 β1 = u+ ,
dove i parametri s± sono delle velocità e legate alla variazione dei bordi tramite
le equazioni
x± = s± t.
(2.156)
Le condizioni (2.155) implicano che nella zona compresa tra x− (t) e x+ (t), la
zona interna, gli invarianti di Riemann β3 e β1 sono costanti e fissati rispettivamente a u− e u+ , mentre, l’unico invariante di Riemann variabile è β2 .
In questa regione, le equazioni di Whitham (2.152) hanno come soluzione
β1 = cost = u+
β3 = cost = u−
s = V2 (u− , β2 , u+ ).
(2.157)
Sostituendo nella (2.135) i valori di u− e u+ si ricava la soluzione per V2
2
(1 − m)m∆
x
1 −
(u + 2u+ + m∆) −
=
3
3 E(m)/K(m)) − (1 − m)
t
(2.158)
dove ∆ = u− − u+ è l’ampiezza della discontinuità iniziale e m = (β2 − u+ )/∆.
La soluzione ottenuta è detta ”Simple-Wave”, soluzione con due invarianti di
Riemann fissati (β1 e β3 )e uno variabile (β2 ).
2
1.8
1.6
1.4
u(x)
1.2
1
0.8
0.6
0.4
0.2
0
−6
−4
−2
0
2
4
6
8
x
Figura 2.6: Dispersive shock wave (DSW) nella KdV per u− = 1 e u+ = 0.
2.4. COSTRUZIONE DI GUREVICH-PITAEVSKII
24
Le velocità auto similari s+ e s− si ottengono ponendo direttamente m = 0
e m = 1 nella (2.158)
s− = s(m = 0) = u+ − ∆,
2
s+ = s(m = 1) = u− + ∆
3
(2.159)
É possibile ricavare il numero d’onda nel limite armonico e l’ampiezza del
solitone nel limite solitonico imponendo le condizioni (2.147) alla (2.50) e (2.55)
.
.
2
−
∆
e
a+ = 2∆.
(2.160)
k =
3
Capitolo 3
Teoria di Whitham in
sistemi non integrabili
La regolarizzazione degli shock per le equazioni non integrabili è un problema aperto e di non semplice soluzione perchè la rottura dell’integrabiltà implica
l’esistenza di un numero finito di leggi di conservazione. Tuttavia, il calcolo
dei parametri s+ , s− e l’ampiezza del solitone as possono essere effettuati dallo
studio delle leggi non mediate al ”trailing edge” e al ”leading edge” in quanto
godono di particolari proprietà. Generalmente due zeri del polinomio F (u, Ai )
collassano semplificando le equazioni di modulazione. Pur non essendo in grado di calcolare gli invarianti di Riemann è possibile riformulare il problema di
Gurevich-Pitaevskii determinando i parametri principali dello shock partendo
dalla forma non diagonale delle equazioni di Whitham.
3.1
Formulazione non integrabile per la KdV
La formulazione della costruzione di Gurevich-Pitaevskii per le equazioni
non integrabili è introdotta a partire dalla KdV di cui sono stati già mostrati i
risultati. Le equazioni di modulazione della KdV sono
%
∂u
∂ $ u2
+
+ uθθ = 0,
∂t
∂x 2
%
u2
∂ $ u2 %
∂ $ u3
+
− θ + uuθθ = 0,
∂t 2
∂x 3
2
(3.1)
(3.2)
∂k ∂ω
+
= 0,
(3.3)
∂t
∂x
dove la mediazione è quella definita nella (2.65) sull’onda di propagazione definita dalla (2.46).
Al ”trailing edge” e al ”leading edge” si può verificare che [7]
F (u)u2 =u3 = F (u3 ),
F (u)u2 =u1 = F (u1 ).
Questo implica automaticamente che
25
(3.4)
3.1. FORMULAZIONE NON INTEGRABILE PER LA KDV
u|u2 =u3 = u3
e
u|u2 =u1 = u1 ,
26
(3.5)
tali relazioni possono essere verificate direttamente dalla (2.70).
Quindi, dalla (3.4) e dalla (3.5) si ricava
F (u)u2 =u3 = F (u),
F (u)u2 =u1 = F (u).
(3.6)
Quest’ultima relazione implica che nei due limiti si ha, u2 = u2 . La prima
equazione di modulazione, nel limite non dispersivo, diventa
ut + uux = 0,
(3.7)
l’equazione di Hopf in u. Nel paragrafo $2.5 è stato esposto il limite semiclassico per cui al ”trailing edge” l’equazione esterna è l’equazione di Hopf, quindi,
la costruzione appena esposta è coerente con i risultati ottenuti nella (2.144) e
(2.146). In definitiva, al ”trailing edge” e al ”leading edge” una equazione di
modulazione collassa sempre nell’equazione di Hopf. Fissato questo concetto si
può formulare la struttura della teoria.
Le equazioni di modulazione (3.1), (3.4) si possono scrivere nella generica
forma differenziale
d
d
Pi (u1 , u2 , u3 ) +
Qi (u1 , u2 , u3 ) = 0.
dt
dx
(3.8)
dove u1 , u2 e u3 sono gli zeri del polinomio G(u). La (3.8) può essere espressa
nella forma
Ut + A(U )Ux = 0
(3.9)
dove U = (u1 , u2 , u3 )T e A(U ) è la matrice dei coefficienti. Si vuole esprimere
la (3.8) in forma caratteristica.
3
5
j=1
(m) dm
cj
dt
U j = 0,
∂
dm
∂
≡
+ Vm ,
dt
∂t
∂x
m = 1, 2, 3,
(3.10)
Ottenuta diagonalizzando A
c(m) A = Vm c(m) ,
(3.11)
dove le Vi , velocità caratteristiche, sono gli autovalori associati al problema.
L’indice m è associato all’equazione caratteristica considerata, mentre j scorre
tra le componente del vettore U e c.
Gli autovettori cj corrispondenti a differenti autovalori Vj sono linearmente indipendenti. Inoltre, se la matrice A è reale e simmetrica tutti gli autovalori sono
reali e ortogonali.
Introducendo un nuovo set di variabili Y = (u, k, a)T , tramite una trasformazione U → Y , si ottiene un nuovo sistema in forma matriciale analogo al
primo
Yt + A(U (Y ))Yx = 0
(3.12)
dove A(U (Y )) è la nuova matrice dei coefficienti. In forma caratteristica diventa
3
5
j=1
(m) dm
bj
dt
Yi =0
m = 1, 2, 3,
(3.13)
3.1. FORMULAZIONE NON INTEGRABILE PER LA KDV
27
dove le velocità caratteristiche Vj (Y ) sono gli autovalori della matrice A(Y ), con
b(m) autovettore applicato da sinistra. Si assuma che autovalori Vm siano reali
e ordinati nel seguente modo: V1 (Y ) < V2 (Y ) < V3 (Y ). In queste condizioni il
sistema definito nella (3.13) diventa iperbolico.
Entrambe le rappresentazioni (3.10) e (3.13) devono possedere le stesse proprietà a cui sono soggette le equazioni di Whitham.
Al ”trailing edge” m = 0, u2 = u3 e a = u3 − u2 = 0, mentre, al ”leading edge”
m = 1, quindi u2 = u1 e k = 0. In entrambi i limiti, nella descrizione con il set
di variabili U , due zeri collassano e il sistema diventa degenere. Un autovalore
della matrice A(U ) acquista doppia molteplicità e due velocità caratteristiche
Vm coincidono.
Un comportamento analogo si avrà anche nella rappresentazione del sistema con
il set di variabili Y , tuttavia nella rappresentazione caratteristica vi sono importanti ripercussioni. Dato che nei due limiti si ha a = 0 o k = 0 un differenziale
scompare automaticamente e il sistema (3.13) ammette una riduzione.
Applicando la costruzione di Gurevich-Pitaevskii, dato che deve essere mantenuto l’ordinamento delle velocità Vi , gli autovalori che collassano devono essere identificati con le velocità dell’onda cnoidale nei limiti armonico e solitonico,
mentre il rimanente autovalore è la velocità u dell’equazione di Hopf.
Se, dopo la riduzione, è nota l’espressione della forma caratteristica si può stabilire una corrispondenza biunivoca tra le restanti componenti di Y . In questi
termini, per la KdV, la velocità dei bordi (2.148) e (2.149) diventano
6
dx−
= V2 (u, k, a = 0)6x=x− ,
dt
6
dx+
= V2 (u, k = 0, a)6x=x+ .
dt
(3.14)
Nel caso di un decadimento da discontinuità iniziale del tipo (2.153), la costruzione di Gurevich-Pitaevskii deve considerare al traling edge si hanno onde
lineari caratterizzate da a = 0 mentre al ”leading edge” si ha il solitone con
k = 0 e i bordi devono raccordarsi con la discontinuità finita
x = s− t :
x = s+ t :
a = 0,
k = 0,
u = u−
u = u+
(3.15)
s+ = V2 (u+ , a+ ),
(3.16)
dove, le velocità s± sono
s− = V2 (u− , k − ),
e a+ , k − sono ottenuti tramite la soluzione del problema di GurevichPitaevskii.
Questo apparto teorico, nonostante non descriva la dinamica dell’onda cnoidale
nella zona interna, può predire esattamente i valori di a+ e k − e, quindi, le
velocità s± al ”trailing edge” e al ”leading edge” senza il calcolo degli invarianti
di Riemann del sistema mediato.
3.1.1
Trailing Edge
Segue dal sistema definito dalla (3.13) che il differenziale lungo la i-esima
caratteristica non è indipendente, ma le componenti del vettore Y sono legate
dalla relazione:
(i)
(i)
(i)
b1 (u, k, a)di u + b2 (u, k, a)di k + b3 (u, k, a)di a = 0.
(3.17)
28
3.1. FORMULAZIONE NON INTEGRABILE PER LA KDV
dove di = d/dt + Vi d/dx. Posto a = 0, un differenziale è nullo e la caratteristica
ammette una riduzione in
(i)
(i)
b1 (u, k, a = 0)di u + b2 (u, k, a = 0)di k = 0,
(3.18)
un sistema 2x2 nelle due variabili (u, k). La risoluzione della (3.18) comporta la
diagonalizzazione del sistema 2x2 e quindi il calcolo degli invarianti di Riemann
del sistema ridotto. Da sottolineare che l’esistenza di questi invarianti non comporta l’esistenza degli invarianti di Riemann nel sistema generale. Le equazioni
diagonalizzate definiscono una relazione locale tra i coefficienti del sistema di
Whitham ridotto.
Se al ”trailing edge” a = u3 − u2 = 0, come precedentemente dimostrato, si
ha
F (u1 , u2 , u3 )|a=0 = F (u3 ) = F (u).
(3.19)
Inoltre, si verifica immediatamente che
lim u2 = u2 ,
(3.20)
lim ω = lim ck = k(c|a=0 ) = k(u3 − k 2 ) = ku − k 3 .
(3.21)
a→0
e che per ω, usando la (2.51),
a→0
a→0
Lo stesso risultato è ottenibile direttamente dalla linearizzazione della KdV
ignorando l’integrazione con le onde cnoidali. Infatti, sostituendo
u = u0 + u1 ei(ωt−kx)
(3.22)
nella (2.23), con |u1 | $| u0 |, si ottiene:
iωu1 − iku0 u1 − iu1 k 3 = 0
e quindi la (2.53).
Il sistema di Whitham ridotto diventa
!
ut + uux = 0,
(k)t + (ω0 (u, k))x = 0.
(3.23)
(3.24)
Si sottolinea che per la KdV il sistema ottenuto è già in forma diagonale
comportando una notevole semplificazione del problema. Procedendo con la
sostituzione k = k(u) si calcola la relazione sulla i-esima caratteristica
dk
dk ∂ω
∂ω
ut +
ux +
ux = 0
du
du ∂k
∂u
(3.25)
sostituendo nella prima delle (3.29) si ottiene una equazione differenziale in
∂ω0
dk
∂u
=
.
0
du
u − ∂ω
∂k
dk
du
(3.26)
Se vale la (2.53), l’equazione diventa
dk
k
1
= 2 =
du
3k
3k
(3.27)
3.1. FORMULAZIONE NON INTEGRABILE PER LA KDV
la cui integrazione da
29
.
2
(u + C1 )
(3.28)
3
L’equazione (3.34) pone la relazione tra i valori di k e u sulla caratteristica
del sistema di Whitham quando l’ampiezza è nulla. La costante C1 è ricavabile
dal raccordo con le condizioni iniziali legate al problema di Gurevich-Pitaevskii.
! −
u = u−
x = s− t
a=0
(3.29)
+
+
x =s t
k=0
u = u+
k=
Dalla seconda condizione si ricava
.
2 +
0=
(u + C1 )
3
dalla prima invece
.
k=
2 −
(u − u+ ) =
3
.
→
2
∆ ≡ k−
3
C1 = −u+
→
k =
−
(3.30)
.
2
∆
3
(3.31)
con ∆ = u− − u+ . In accordo con la (3.16) si deve calcolare la velocità auto
similare s− . Le velocità caratteristiche del sistema di Whitham ridotto derivano
dal confronto con le equazioni (3.24) e (3.25), quindi,
∂ω0
= u − 3k 2 .
∂k
u,
(3.32)
La seconda espressione indica la velocità di gruppo di un pacchetto di onde derivante dalla legge di relazione di dispersione (3.21). Per una completo
accordo tra la teoria di Whitham e la formulazione non integrabile si devono identificare le velocità caratteristiche (3.32) con le velocità dal problema di
Gurevich-Pitaevskii integrabile (2.143). Dato che la KdV è una equazione completamente integrabile tale identificazione è abbastanza agevole.
Partendo da u,
E(m)
u = u3 + 2(u3 − u1 )
,
(3.33)
K(m)
tramite la definizione di k, nella (2.46), si ricava
√ √
2 3 u3 − u1
K(m) =
.
π
k
(3.34)
Sostituendo nella (3.66) si ottiene
u = u1 +
Per a = 0, u = u3 ,
ricavando
√
(u − u1 ) =
√
u3 − u1 3k.
(3.35)
√
u − u1 3k,
(3.36)
√
u1 = u − 3k 2 .
(3.37)
V2 (u, k, a = 0) = ∂ω0 /∂k
(3.38)
Questa espressione coincide con la velocità di gruppo nel limite lineare dell’onda
cnoidale, ricavabile dalla (2.53). Con questa condizione la velocità caratteristica
V2 può essere identificata con
3.1. FORMULAZIONE NON INTEGRABILE PER LA KDV
30
e s−
∂ω0 (u− , k − )
= u+ − ∆.
(3.39)
∂k
Risultato in accordo con quello calcolato nella costruzione di GurevichPitaevskii (2.159).
s− =
3.1.2
Leading Edge
Allo stesso modo, per il ”leading edge”, si procede integrazione indiretta
della i-esima caratteristica con k = 0.
(i)
(i)
b1 (u, 0, a)di u + b2 (u, 0, a)di a = 0
(3.40)
Cercando la relazione locale a = a(u) si calcola, come per il ”trailing edge”, l’integrale della (3.40) e quindi s+ . Sulla base degli elementi esposti nel
paragrafo $ 3.1 introducendo un opportuno set di variabili si può semplificare il
problema, in particolare si introduce il numero d’onda coniugato
$"
u2
du #
1
%−1
1
π (u3 − u1 ) 2
= √
,
(3.41)
2 3 K(m" )
G(u)
u1
#
con m" = 1 − m. Si nota che se m" = 0, k̃ = as /3. Si definisce, inoltre, la
quantità
k
Λ= .
(3.42)
k̃
k̃ = π
Il vettore Y (u, k, a) può essere sostituito dal nuovo vettore Y " = Y " (u, k̃, Λ). Le
nuove quantità k̃ e Λ possiedono stesse le proprietà, rispettivamente, di a e k.
m = 1 → m" = 0 →
m = 0 → m" = 1 →
Λ =0
k̃ = 0
(3.43)
Con queste nuove variabili le condizioni di Gurevich-Pitaevskii diventano
! −
x = s− t
k̃ = 0
u = u−
(3.44)
+
+
x =s t
Λ =0
u = u+
Per le equazioni di modulazione nel limite solitonico si ha che la prima decade
nell’equazione di Hopf, mentre per la seconda equazione si ha che
ω → ω̃ = ω̃(u, k̃, λ) = ck̃.
Usando (3.42) e la (3.45) si ricava
$ %
$ %
$
%
k̃Λ + ω̃ = 0 → kΛt + ω̃Λx + Λ kt + ω˜x = 0.
t
x
(3.45)
(3.46)
Equivalentemente al ”trailing edge”, se λ = 0 su una caratteristica i valori di k̃
e u devono essere connessi da una relazione del tipo k̃ = k̃(u) che non dipende
dalla soluzione particolare. Considerando l’equazione (3.46) nell’intorno della
caratteristica menzionata quando k $ 1, per il termine tra parentesi nella (3.45)
si ha
%
dk̃ ∂u ∂ ω˜s dk̃
dk̃ $ ∂ ω˜s
+
ux = ux
−u .
(3.47)
dũ ∂t
∂k du
du ∂u
3.1. FORMULAZIONE NON INTEGRABILE PER LA KDV
31
dove indichiamo con ω̃s = ω̃(u, 0, k̃) la relazione che può essere chiamata relazione di dispersione solitonica.. Quindi, sostituendo nella (3.45),
λt +
% ∂ ω̃ ω˜s
λ , dk̃ $ ∂ ω̃s
∂λ
s
+ ux
−u +
= O(λ )
du
∂u
∂u
∂x
k̃
k̃
(3.48)
La caratteristica da la relazione k̃ = k̃(u) e non deve dipendere da λ, quindi,
l’espressione tra parentesi quadre deve essere identicamente nulla. Si ricava, in
definitiva,
∂ ω̃s
dk̃
∂u
=
(3.49)
du
u − ∂ ω̃0
∂ k̃
e
ω˜s
∂λ
λx = O(λ ).
(3.50)
∂x
k̃
Notiamo che le espressioni (3.50) e (3.26) sono identiche a meno di sostituire
l’opportuna legge di dispersione. Le velocità caratteristiche in questo limite
corrispondono a
ω̃s
,
(3.51)
u,
k̃
dove il primo termine è la velocità legata alla equazione di Hopf, mentre il secondo è la velocità associata alla doppia caratteristica ottenuta dalla riduzione.
Integrando la (3.49) con la legge di dispersione solitonica (3.45) e imponendo le
condizioni (3.44) si può ottenere il valore del parametro s+ . La costruzione con
il set di variabili Y " conduce al calcolo di s+ e quindi alla soluzione del problema di Gurevich-Pitaevskii. Si osserva, innanzitutto, che k̃ e k possono essere
interpretati come numero d’onda e frequenza della onda di propagazione coniugata associata alle stesse radici uj , ma con curva di potenziale G(u) invertita.
Questi risultati si ottengono automaticamente definendo una trasformazione di
coordinate del tipo
x → ix̃
t → it̃
u → ũ.
(3.52)
λt +
In questo modo la KdV si trasforma in una nuova equazione che porta segno
opposto nel termine dispersivo
ũt̃ − ũx̃x̃x̃ + ũũx̃ .
(3.53)
Cercando una soluzione di propagazione ũ = ũ(θ̃) = ũ(x̃ − ct̃), chiamata onda
coniugata a cui è associata la stessa velocità di propagazione c del caso reale, si
ottiene agevolmente
(ũ)2θ̃ = G(ũ),
con ũ(θ + L̃) = ũ(θ)
(3.54)
Mentre k̃ = 2π/L è definito nella (3.42)
Nel caso particolare della KdV le soluzioni di propagazione sono in un caso
u = u(θ) e nell’altro iũ = (−iθ̃). Entrambe hanno lo stesso set di radici uj
e rappresentano proiezioni sull’asse reale e immaginario della stessa funzione
ellittica. La differenza fondamentale sta che per le due onde il numero d’onda è
simile, ma con peculiarità diverse.
Infatti, per l’onda coniugata il limite armonico cade per u1 = u2 , mentre
il limite solitonico per u3 = u2 , viceversa per il caso reale. Quindi il limite
3.1. FORMULAZIONE NON INTEGRABILE PER LA KDV
32
armonico dell’onda coniugata corrisponde al limite solitonico dell’onda reale e
per entrambe le equazioni il limite armonico corrisponde alla propagazione di
onde lineari con una specifica legge di dispersione. Senza la richiesta di integrabilità, la legge di dispersione per la KdV reale e coniugata si trova mediante
linearizzazione attorno ad una valore fisso di fondo.
Nello specifico, la KdV coniugata si linearizza, analogamente alla (3.22),
attorno ad un valore u0 = +u, ottenendo
ω˜s = k̃+u, + k̃ 3 .
(3.55)
Equivalentemente, si definisce la media sulla scala veloce, della famiglia di
soluzioni definite nella (3.53). Per un osservabile F si ha
+F ,(u1 , u2 , u3 ) =
k̃
π
"
u2
u1
F (ũ)
dũ #
,
G(ũ)
dove +, è la media. Con l’integrale definito nella (3.56), si mostra che
6
6
+F ,(u1 , u2 , u3 )6
= F (u1 ),
u2 =u1
quindi, tramite la (3.5),
+F (u),k=0 = F (u)k=0 ,
(3.56)
(3.57)
(3.58)
in particolare si ha che +u, = u. Questa proprietà è legata alla forma quadratica
della G(u) nell’intorno del limite armonico e solitonico.
Segue dalla (3.58) che
ω˜s (+u,, k̃) = ω˜s (u, k̃).
(3.59)
Usando la (3.55), si ottiene per la KdV coniugata, quando k = 0,
ω˜s (u, k̃) = k̃u + k̃ 3 .
(3.60)
Il problema ora è formalmente risolto perché tramite imposizione delle (3.44)
si ottiene con passaggi identici al caso reale, ma con la legge di dispersione
leggermente diversa
.
2
k̃ =
(C2 − u),
(3.61)
3
dove C2 costante di integrazione. Va sottolineato che per il caso reale, tramite
il cambiamento formale k → ik̃ e ω0 → iω˜0 , la soluzione del caso solitonico la
si può ottenere direttamente dalla legge dispersione lineare. La soluzione (3.61)
è una relazione generale tra k̃ e u sulla caratteristica del sistema di Whitham
quando λ = 0 e non è una soluzione globale. Ora, tramite la costruzione di
Gurevich-Pitaevskii si ricava il valore delle costanti. Imponendo la (3.44) si
ottiene
.
2
0=
(C2 − u− ) → C2 = u− .
(3.62)
3
Imponendo la seconda otteniamo il valore del numero d’onda coniugato:
.
2
+
k̃ =
∆.
(3.63)
3
3.2. FORMULAZIONE NON INTEGRABILE PER SISTEMI
BIDIREZIONALI
33
Si osserva, dalla (3.31), che k̃ + = k − . Seguendo ancora la costruzione del
”trailing edge” della KdV reale, la velocità s+ si ottiene dalla velocità di gruppo
associata alla relazione di dispersione coniugata (3.60)
s+ =
ωs (u+ , k̃ + )
2
= u+ + ∆.
3
k̃
(3.64)
In accordo con la soluzione globale del sistema nel limite solitonico (2.159).
L’ampiezza as del solitone si ottiene dalla relazione (2.58) che per l’onda coniugata diventa
2
s+ = u+ + a+ ,
(3.65)
3
usando la (3.64) si ricava
a+ = 2∆ = 2(u− − u+ ).
(3.66)
in accordo con la (2.160).
3.2
Formulazione non integrabile per sistemi bidirezionali
La formulazione di Gurevich-Pitaevskii si può estendere ai sistemi bidirezionali che vantano numerose applicazioni in ottica e in meccanica dei fluidi. Più
precisamente a sistemi 2x2 strettamente iperbolici modificati da una piccola dispersione. Si consideri un generico sistema bidirezionale espresso nella simbolica
forma:
φt = K4 (φ, φx , φxx , ...),
(3.67)
dove φ è un vettore a due componenti e K4 è un operatore differenziale che
mescola le derivare spaziali φx . Tale che la (3.67) ammetta una legge di dispersione lineare ω = ω0 (k) che nel limite non dispersivo diventi ω → k se
k $ 1. Si assume che il sistema (3.67) ammetta almeno tre leggi di conservazione
indipendenti
∂Qi
∂Pi
+
= 0,
i = 1, 2, 3.
(3.68)
∂t
∂x
dove Pi sono le quantità conservate e Qi i corrispondenti flussi. In aggiunta alle
equazioni (3.68) si deve considerare anche la conservazione del numero d’onda.
∂k ∂ω
+
=0
∂t
∂x
(3.69)
Generalmente i sistemi fisici considerati, meccanica dei fluidi e ottica, conservano densità ρ e momento ρu. In questa breve esposizione si tratterà un modello che descrive la dinamica di un gas dove due leggi di conservazione P1,2
corrispondano, nel limite non dispersivo , alle equazioni di Eulero
(ρ)t + (ρu)x = 0,
(ρu)t + (ρu2 + p(ρ))x = 0.
(3.70)
dove p(ρ) è la pressione del gas. Assunzione generale per tali modelli è che
il limite dispersivo ammetta la rappresentazione in funzione degli invarianti di
3.2. FORMULAZIONE NON INTEGRABILE PER SISTEMI
BIDIREZIONALI
34
Riemann. Si cerca per il sistema (3.67) una soluzione di propagazione, periodica
e a singola fase
φ(x, t) = φ(θ)
con
φ(θ + 2π/k) = φ(θ),
(3.71)
dove θ = x − ct. La soluzione deve essere espressa in funzione delle variabili ρ,
u, k e a. Dove ρ e u la densità e la velocità media del fluido, mentre k e a sono
gli usuali numero d’onda e ampiezza della soluzione periodica. Come ultima
assunzione si ammette che la riduzione del tipo (2.14) conduca ad un potenziale
Q(ρ) che possieda le proprietà discusse nella sezione precedente, ovvero una
riduzione lineare se a = 0 e iperbolica se k = 0. I tipici sistemi a possedere
queste caratteristiche sono la NLS [10] o il sistema Kaup-Boussinesq [11]. Le
tipiche condizioni iniziali a gradino per un sistema del tipo (3.70) sono
!
!
ρ = ρ−
ρ = ρ−
x>0
(3.72)
t=0:
− x < 0,
u = u−
u=u
dove ρ± e u = u± sono costanti. Studi analitici sui decadimenti da discontinuità iniziale in sistemi dispersivi integrabili mostrano che la soluzione asintotica
consiste in di tre regioni separate da due onde in propagazione: un’onda di rarefazione e/o onda di shock. La descrizione dello shock dispersivo generalmente si
esegue tramite la teoria di Whitham, il sui sistema si può esprimere nella forma:
∂
∂
P j (ρ̄, ū, k, a) +
Q (ρ̄, ū, k, a) = 0 ,
∂t
∂x j
j = 1, 2, 3 ,
(3.73)
∂
∂
k+
ω(ρ̄, ū, k, a) = 0 ,
(3.74)
∂t
∂x
dove ρ̄ and ū sono la densità e la velocità mediate sulla famiglia si soluzioni
(3.70), k è il numero d’onda e a è l’ampiezza. Usando gli argomenti esposti
nella sezione precedente, il sistema di Whitham ammette una riduzione esatta
per a = 0 e k = 0 che deve corrispondere con il limite non dispersivo del sistema
(equazioni di Eulero mediate)
7 8
7 8
ρ t + ρ u x = 0,
7 8
7
8
(3.75)
ρu t + ρ u2 + p(ρ) x = 0,
Ovviamente il sistema (3.75) può essere ottenuto direttamente dal sistema di
Whitham considerando i limiti lineare e iperbolico. La legge di conservazione
del numero d’onda (3.69), nel limite lineare assume la forma
a = 0,
kt + (ω0 (ρ̄, ū, k))x = 0
(3.76)
dove ω = ω0 (ρ̄, ū, k) è la legge di dispersione lineare del sistema ottenuta considerando piccole perturbazioni su un fondo ρ = ρ̄, u = ū.
Le condizioni iniziali per i sistemi bidirezionali nel caso di decadimento da una
discontinuità iniziali si ottengono per analogia con le condizioni al contorno
(3.15) della KdV
x = s− t :
x = s+ t :
a = 0 , ρ̄ = ρ− , ū = u− ,
k = 0 , ρ̄ = ρ+ , ū = u+ ,
(3.77)
3.2. FORMULAZIONE NON INTEGRABILE PER SISTEMI
BIDIREZIONALI
35
dove la dipendenza delle velocità dei bordi s± in funzione dei parametri iniziali
deve essere trovata. La determinazione delle velocità si ottiene dalla richiesta di
una relazione locale tra k e u lungo le curve caratteristiche definite dalla velocità del sistema di Whitham ridotto (a = 0 e k = 0). Nei sistemi bidirezionali,
tipicamente, si deve imporre la relazione tra la soluzione del modello non dispersivo (3.75) e l’espressione della legge di conservazione del numero d’onda (3.74)
espressa in funzione delle variabili mediate. In definitiva si cerva la relazione
tra le componenti ρ e u in funzione dei parametri iniziali
Φ(ρ+ , ρ− , u+ , u− ) = 0
(3.78)
Da sottolineare che negli shock conservativi, questo metodo è applicabile finché
la pressione p è funzione della sola densità ρ.
Capitolo 4
Nonlinear Schrödinger
quintica
In questo capitolo si vuole applicare la costruzione di Gurevich-Pitaevskii
per la nonlinear Schrödinger con un termine quintico perturbativo (CQNLS) in
regime semiclassino
i#ψz +
#2
α
ψxx − |ψ|2 ψ + |ψ|4 ψ = 0,
2
2
(4.1)
dove ψ rappresenta un campo complesso, 0 < # $ 1 parametro per la dispersione e il termine α è un parametro esterno di controllo per la nonlinearità quintica.
L’equazione (4.1) può essere considerata la generalizzazione della nonlinear
Schrödinger (NLS) iψz + ψxx + a|ψ|2 ψ = 0. Per a < 0 la NLS descrive un
modello con nonlinearità ”defocusing”, viceversa, se a > 0 si parla di caso
”focusing”.
Quando a < 0 la soluzione che descrive del modello è del tipo
#
ψ =
−2/a sech(x) exp (iz). Come la NLS, la CQNLS è un’equazione universale che descrive situazioni fisiche diverse. I modelli più importanti a cui è
associata la CQNLS sono lo studio dei plasmi [12], ottica non lineare [13] e condensazione di Bose-Einstein [14]. Nella condensazione di Bose, le non linearità
cubica e quintica rappresentano rispettivamente l’interazione a due e tre corpi.
In ottica, la CQNLS descrive la propagazione di un’onda elettromagnetica stazionaria in un materiale con una non linearità Kerr di tipo saturabile troncata al
primo ordine. Questa classe di equazioni può essere ulteriormente generalizzata
introducendo una funzione g(|ψ|2 ), per la parte non lineare, che vari a seconda
del modello considerato. Per la CQNLS è g(|ψ|2 ) = |ψ|2 − α2 |ψ|4 . L’equazione
associata alla g(ψ) è la nonlinear Schrödinger generalizzata (GNLS)
i#ψz +
#2
ψxx − g(|ψ|2 )ψ = 0.
2
(4.2)
Spesso è conveniente risolvere i modelli in funzione di g(|ψ|2 ) per trovare caratteristiche più generali della soluzione. A differenza della NLS, la CQNLS non
appartiene alla classe di equazioni integrabili e quindi non può essere risolta mediante il metodo della trasformata spettrale [4]. Tuttavia, la CQNLS ammette
soluzioni per onda solitaria, legate da un meccanismo di bistabilità. Inoltre, la
CQNLS è una equazione che presenta fenomeni di shock, e, data la sua natura
36
4.1. TEORIA ANALITICA DELLO SHOCK GENERATO DA UNA
DISCONTINUITÀ INIZIALE
37
non integrabile la risoluzione non può avvenire tramite la teoria di Whitham,
bensı̀, tramite la riformulazione del problema di Gurevich-Pitaevskii per modelli
non integrabili presentata nel precedente capitolo.
In linea con la trattazione generale degli shock, comunemente usata per la NLS,
si procede applicando la trasformazione di Madelung
* "
+
#
i
ψ = ρ(x, z) exp
u(x, z)dx ,
(4.3)
#
dove ρ > 0 rappresenta l’intensità del campo ψ. Sostituendo nella (4.2) si ottiene
l’equivalente idrodinamico del modello
ρz + (uρ)x = 0,
$ρ
ρ2 %
xx
uz + uux + f " (ρ)ρx = #2
− x2 ,
4ρ
8ρ x
(4.4a)
(4.4b)
dove f (ρ) è la non linearità della CQNLS
f (ρ) = ρ −
α 2
ρ .
2
(4.5)
Nell’equivalente idrodinamico, ρ può essere interpretata come la densità del fluido e la (4.5) è il corrispondente delle equazioni di stato. La funzione u(r, z),
il valore locale del vettore d’onda, è la velocità del fluido. La prima equazione
delle (4.4) rappresenta la conservazione della massa e la seconda la conservazione dell’impulso e i termini al secondo membro della (4.4b) sono associati alla
dispersione. Come per la KdV, si dovrà studiare il sistema al variare di # e la
dispersione avrà un ruolo determinante solo per tempi confrontabili con il tempo
di rottura tc . Il caso # = 0 corrisponde al limite non dispersivo della CQNLS
ρz + (uρ)x = 0
(4.6a)
uz + uux + f " (ρ)ρx = 0
(4.6b)
Si sottolinea che il limite α → 0, corrispondente alla NLS e deve essere valido
ad ogni livello dell’analisi dello shock.
4.1
Teoria analitica dello shock generato da una
discontinuità iniziale
Si consideri il problema descritto dalla (4.4), con f (ρ) definita nella (4.5),
relativo alle condizioni iniziali per ρ e u del tipo:
!
ρ0 x < 0
ρ=
u(x, 0) = 0,
(4.7)
1 x>0
dove ρ0 > 1. Si assume che il valore iniziale della u sia nullo su tutto lo spazio.
All’inizio dell’evoluzione, si generano onde di propagazione lineari in accordo con
la legge di dispersione ottenuta tramite la linearizzazione di (4.4) nell’intorno
di ρ = ρ0 e u = u0 , valori di fondo del pacchetto d’onda e # ∼ 1. Sostituendo:
ρ = ρ0 + ρ1 exp [i(kx − ωz)] e u = u0 + u1 exp [i(kx − ωz)] con |ρ1 | $ 1, si
ottiene, trascurando i termini di ordine superiore, l’espressione della legge di
dispersione. Dalla (4.4a), sostituendo l’espressione di ρ e u si ricava:
4.1. TEORIA ANALITICA DELLO SHOCK GENERATO DA UNA
DISCONTINUITÀ INIZIALE
(ω − ku0 ) = k
ρ0 u1
,
ρ1
38
(4.8)
mentre, dalla (4.4b), trascurando i termini ρ21 , si ottiene
(ω − ku0 ) + f " (ρ0 )k
ρ1
k 3 ρ1
=−
.
u1
4 u1
(4.9)
Sostituendo la (4.9) nella (4.8) si ricava la legge di dispersione ω = ω(k):
.
k2
ω = ω(k, u0 , ρ0 ) = ku0 ± k f " (ρ0 )ρ0 + .
(4.10)
4
Per la f (ρ) data dalla (4.5) si ricava la legge di dispersione per la CQNLS
.
k2
ω = ω(k, u0 , ρ0 ) = ku0 ± k ρ0 − αρ20 +
(4.11)
4
4.1.1
Limite non dispersivo
Si vuole calcolare la soluzione del sistema (4.6), per le condizioni iniziali (4.9).
Data la bidirezionalità del sistema, generalmente il problema si risolve con una
combinazione di due onde in propagazione in direzioni opposte. Introducendo
gli invarianti di Riemann (vedi Appendice A)
λ± = u ± R(ρ),
con
#
1
R(ρ) = ρ − αρ2 − √ arctan
α
Il sistema (4.6) diventa
(4.12)
.
1
− 1,
αρ
∂λ±
∂λ±
+ V±
=0
∂z
∂x
dove V±
V± = u ±
#
ρ − αρ2 ,
(4.13)
(4.14)
(4.15)
sono le velocità caratteristiche. Da sottolineare che la condizione di sistema
iperbolico è valida finché le velocità sono reali, questo è vero se α < 1/ρ0 .
Nel limite α → 0, si devono ottenere gli invarianti di Riemann della NLS.
Sviluppando la (4.12) e la (4.15) per α piccoli, si ottengono le espressioni degli
invarianti di Riemann e delle velocità caratteristiche per la NLS [10]:
√
r± = u ± 2 ρ,
V± =
3
1
r± + r∓ .
4
4
(4.16)
Dato che il sistema è invariante per trasformazioni del tipo x → Cx e z → Cz,
le (4.14) si possono esprimere in funzione della variabile autosimilare ζ = x/z
$
%
d
V+ − ζ dζ
λ+ = 0
$
%
(4.17)
d
V− − ζ dζ
λ− = 0.
4.1. TEORIA ANALITICA DELLO SHOCK GENERATO DA UNA
DISCONTINUITÀ INIZIALE
39
Si cercano soluzioni di tipo ”simple-wave” caratterizzate da un invariante di
Riemann costante e l’altro variabile. Le possibili soluzioni del sistema (4.17)
sono
V+ = ζ
λ− = λ0− = cost,
(4.18)
oppure
V− = ζ
λ+ = λ0+ = cost.
(4.19)
Dato che nelle condizioni iniziali si ha che ρ = ρ0 e u = 0, per x < 0, si sceglie
come soluzione la (4.19). Imponendo le condizioni (4.7) si ha per λ0+ :
λ0+ =
9
1
ρ0 − αρ20 − √ , arctan
α
.
1
− 1 = R(ρ0 ),
αρ0
(4.20)
che, espresso in funzione di u diventa
u = R(ρ0 ) − R(ρ).
(4.21)
Sostituendo la (4.21) nella (4.19), si ottiene la soluzione cercata, detta left simple
wave (LSW).
.
.
9
% x
#
1 $
1
1
ρ0 − αρ20 − 2 ρ − αρ2 − √
arctan
− 1 − arctan
−1 = ,
αρ0
αρ
z
α
(4.22)
questa equazione determina l’andamento, in forma implicita, dell’intensità ρ in
funzione della variabile auto-similare ζ = x/z. Per un completo raccordo tra
le condizioni iniziali (4.7) e la soluzione (4.22) si individua un punto x = x−
1
−
tale che: per x < x−
1 si ha ρ = ρ0 = cost, mentre, per x > x1 si ha l’onda di
rarefazione, la soluzione del sistema per x < 0.
Quindi, x = x−
1 è un punto di discontinuità debole in cui il sistema si propaga
con la velocità vs , ottenuta sostituendo ρ = ρ0 nella (4.22), pari alla velocità
9
x−
1
= −vs = − ρ0 − αρ20 .
z
(4.23)
Successivamente, per x > 0 si ha: ρ = 1 e u = 0. Queste condizioni però, non
sono in accordo con l’Eq. (4.21) utilizzata per la costruzione della LSW.
L’onda di che soddisfa le condizioni (4.7) per x > 0, si ottiene in modo analogo
alla (4.22) ed ha la forma:
.
√
1
1
0
√
λ+ = 1 − α −
arctan
− 1 = R(1),
(4.24)
α
α
.
.
% x
#
√
1 $
1
1
2
1 − α − 2 ρ − αρ − √
arctan
− 1 − arctan
− 1 = . (4.25)
α
αρ
z
α
Per connettere le due ”simple-wave”, entrambe soluzioni valide del sistema, si
introduce una distribuzione intermedia
ρ(x/z) = ρ− = cost
u(x/z) = u− = cost
(4.26)
tale da raccordare le due onde di rarefazione. Il raccordo tra la distribuzione
intermedia e la LSW iniziale (4.22) avviene in un secondo punto di discontinuità
4.1. TEORIA ANALITICA DELLO SHOCK GENERATO DA UNA
DISCONTINUITÀ INIZIALE
40
debole x = x+
1 . La connessione tra la distribuzione intermedia e la LSW finale
(4.25) può avvenire esclusivamente tramite una right simple wave (RSW). La
RSW si muove in direzione opposta rispetto alle due soluzioni e ciò comporta
la comparsa di una zona compresa tra due punti di discontinuità debole, x−
2 e
x+
in
cui
il
sistema
assume
valore
multiplo
in
virtù
della
sovrapposizione
della
2
RSW con la distribuzione media (4.26) e la LSW (4.25) valida per x > 0.
La RSW si ottiene imponendo per la (4.18) le condizioni ρ = 1 e u = 0 in modo
da connettere la LSW dove λ0+ = R(1):
.
√
1
1
(4.27)
λ0− = √ arctan
− 1 − 1 − α,
α
α
u = R(ρ) − R(1)
(4.28)
e l’andamento in funzione di ζ
#
√
1 $
2 ρ − αρ2 − 1 − α + √
arctan
α
.
1
− 1 − arctan
αρ
.
% x
1
− 1 = . (4.29)
α
z
Le soluzioni (4.22) e (4.29) devono essere connesse alla distribuzione intermedia ρ = ρ− e u = u− . Imponendo che le velocità u della LSW iniziale e la
RSW siano uguali a u− in un punto di discontinuità x = x+
2 si ricava il valore
del parametro ρ− relativo alla discontinuità iniziale.
! −
u = R(ρ0 ) − R(ρ− )
(4.30)
u− = R(ρ− ) − R(ρ1 )
quindi
2R(ρ− ) = R(ρ0 ) + R(1),
(4.31)
che definisce completamente ρ− . Noto il parametro ρ− , si ricava dalla (4.28) u−
%
1$
u− = R(ρ− ) − R(1) =
R(ρ0 ) − R(1) .
(4.32)
2
−
I punti di discontinuità interni x = x+
1 e x = x2 si ottengono sostituendo i
−
−
valori di ρ e u nelle (4.29) e (4.22),
x+
1
= u− − vs (ρ− ),
z
x−
2
= u− + vs (ρ− )
z
(4.33)
Questi punti corrispondono alle discontinuità deboli con cui l’onda si propaga
con la velocità vs , ma nella direzione opposta. L’intera distribuzione di intensità
è mostrata in figura (4.1a). La soluzione formale del sistema (4.4) porta nella
+
regione x−
2 < x < x2 ad una distribuzione a tre valori e quindi non fisica. Questa
regione indica che si genera uno shock dispersivo nella regione di transizione tra
ρ = ρ− , u = u− e ρ+ = 1, u+ = 0. La struttura che si forma è presentata
qualitativamente in figura (4.1b). Si vuole sottolineare che i bordi della zona
di oscillazione x±
2 non devono necessariamente coincidere nei due casi. La zona
di oscillazione sorge, come nella KdV, come effetto di compensazione tra la
dispersione e gli effetti non lineari. Nel punto x+
2 , il ”leading edge”, l’onda
si trasforma in un pacchetto d’onda lineare, mentre, nel punto x−
2 , il ”trailing
edge”, si ha un solitone.
4.2. SOLITONI DARK E ANTIDARK
41
Figura 4.1: Soluzione del problema di Riemann associato associato alla
CQNLS. (a) soluzione del limite non dispersivo. (b) figura schematica per la
rappresentazione dello shock.
Si può osservare la zona di oscillazione dello shock caratterizzata da due
differenti scale spaziali: l’oscillazione veloce dell’intensità e la variazione lenta
del profilo di modulazione. Questo suggerisce che lo shock può essere trattato
con la teoria di Whitham.
4.2
Solitoni Dark e Antidark
Prima di passare alla descrizione dello shock si devono calcolare le soluzione per onde di traslazione della CQNLS [15] quando i termini dispersivi sono
confrontabili con la non linearità (# = 1). Partendo dalla (4.2), si impone
inizialmente che il campo ψ sia stazionario in t
ψ = u(x, z)eiβz ,
(4.34)
42
4.2. SOLITONI DARK E ANTIDARK
dove β = g(ρb ) con ρb = |ψ0 |2 intensità di fondo. Conseguentemente, anche
|u(x, t)|2 → ρb quando x → ±∞. La GNLS (4.2) espressa in termini di u(x.t)
diventa:
:
;
1
iuz + uxx − g(|u|2 ) − g(ρb ) = 0.
(4.35)
2
La soluzione di solitone si trova imponendo un onda di propagazione, localizzata,
nella forma
u(x, z) = f (θ)eiφ(θ) ,
(4.36)
con f e φ funzioni reali e θ = x − vz. Queste soluzioni dipendono da due
parametri interni la velocità v e l’intensità di fondo ρb . Sostituendo la (4.36)
nella (4.35) si ottengono le condizioni che devono soddisfare f e φ. Con questa
sostituzione le derivate diventano:
∂
∂
=
,
∂x
∂θ
∂
∂
= −v .
∂z
∂θ
Analizzando singolarmente le derivate del campo u si ha
*
+
fθ
iuz = v φθ − i
u,
f
(4.37)
(4.38)
e
1
1
uxx =
2
2
/*
fθ
+iφθ
f
+
+
x
*
fθ
+iφθ
f
+2 0
*
+
fθθ φθθ φ2θ φθ fθ
u=
+i
− +i
u. (4.39)
2f
2
2
f
Sostituendo nella (4.35) e separando parte reale e immaginaria si ha:
f 2 φθθ
+ f fθ φθ ,
2
7
8
fθθ
φθ 8 7
+ vφθ 1 −
− g(ρb ) + g(f 2 ) = 0.
2f
2
vf fθ =
(4.40a)
(4.40b)
Integrando una volta rispetto a θ la (4.40a) si ha
vf 2 = f 2 φθ + 2B,
(4.41)
con B costante. Imponendo la condizione che f 2 → ρb per x → ±∞, dove la
fase φ è stazionaria, si ricava B = ρb /2v. Si ha, infine, per la parte immaginaria
*
+
ρb
φθ = v 1 − 2 .
(4.42)
f
La (4.42) impone una condizione sulla fase φ che sostituita nella (4.40b) determina l’equazione che deve soddisfare f (θ)
*
+ /
0
fθθ
v2
ρ2b
2
+
f − 3 − g(f ) − g(ρb ) f = 0.
(4.43)
2
2
f
Si riconosce nella (4.43) un sistema Hamiltoniano. Sostituendo g(f ) = f 2 −
(α/2)f 4 si ricava
fθθ = −v 2 f +
v 2 ρ2b
∂V
+ 2g(ρb )f + 2f 3 + γf 5 = −
f3
∂f
(4.44)
43
4.2. SOLITONI DARK E ANTIDARK
che definisce il potenziale V (f ) dell’Hamiltoniana
*
+ 2
α 6 1 4
f
v 2 ρ2b
2
V (f ) = f − f + v + 2g(ρb )
.
+
6
2
2
2f 2
(4.45)
In conclusione, l’onda di traslazione deve soddisfare
Ĥ(f, f˙) =
fθ2
+ V (f ).
2
(4.46)
Le soluzioni di solitone corrispondono alle orbite omocline/eterocline della (4.46).
Il solitone dark con intensità di fondo ρb si ricava dalle orbite Hamiltoniane ca√
ratterizzate dal livello E0 = Ĥ(f = ρb , f˙ = 0), tale che E0 − V (f ) ha un
√
doppia radice per f = ρb . Definendo 2k0 ≡ v 2 + 2g(ρb ) e C0 = vρb si ottiene
il livello energetico del solitone dark
α
ρ2
E0 = − ρ3b + b + v 2 ρb .
3
2
(4.47)
Tuttavia, per analisi analitiche successive è più conveniente esprimere il sistema
in funzione dell’intensità ρ ≡ f 2 , che, sostituita nella (4.46) da
4
ρ2θ = − αρ4 + 4ρ3 − 8k0 ρ2 + 8E0 ρ − 4C02 ≡ Q(ρ),
3
(4.48)
Si riconosce un nuovo sistema Hamiltoniano in cui ρ è la posizione e p = ρ̇ ≡
dρ/dθ è il momento, che evolvono secondo le equazioni canoniche
p2
+ V (ρ),
2
*
+
ρ2θ
α
C2
= Ĥ − V (ρ) = 2ρ − ρ3 + ρ2 − 2k0 ρ + 2E0 − 0 ,
2
3
ρ
ṗ = −
∂ Ĥ
,
∂ρ
ρ̇ =
∂ Ĥ
;
∂p
Ĥ(ρ, ρ̇) =
che unita con la (4.47) porta alla forma del potenziale V (ρ)
/
0
C02
α 3
3
2
2
V = 2ρ
(ρ − ρb ) − (ρ − ρb ) + 2k0 (ρ − ρb ) −
.
3
ρ0
(4.49)
(4.50)
(4.51)
I solitoni corrispondono alle orbite omocline della (4.50) che si formano dal
punto di sella dell’Hamiltoniana (p, ρ) = (0, ρb ) e caratterizzate dal livello energetico E #
= V (ρb ) = −2C02 . L’equazione (4.48) deve soddisfare la condizione ρθ =
Q(ρ), quindi sono ammesse soluzioni reali solo nelle zone in cui
Q(ρ) = [E − V (ρ)] > 0. Per convenienza si riporta l’analisi in funzione di
−Q(ρ) in modo da studiare il comportamento dei minimi di V (ρ) − E. Per
α -= 0 il potenziale V (ρ) − E è quartico con una doppia ”buca” a cui è associato
un doppio loop nello spazio delle fasi (ρθ , ρ) come mostrato in figura 4.2. Lo
spazio delle fasi mostra l’esistenza due stati legati caratterizzati rispettivamente
da una soluzione che varia tra per ρm < ρ < ρb , dove ρm è l’intensità di minimo
e una per ρb < ρ < ρa con ρa intensità di picco. Si ha in generale la coesistenza
di due soluzioni corrispondenti ad una coppia di solitone dark e antidark (solitone bright su un piedistallo ρb ). Le intensità ρm < ρb < ρa sono gli zeri reali di
V (ρ) − E e la loro forma esplicita si ricava direttamente dalla (4.48). In questo
limite, limite solitonico, il funzionale Q(ρ) = Ĥ − V (ρ) presenta un doppio zero
44
4.2. SOLITONI DARK E ANTIDARK
2
5
(b)
4
1.5
3
2
1
ρθ
V(ρ)−E
1
0.5
ρm
ρa
ρ
b
0
−1
0
−2
−3
(a)
−0.5
−4
−1
0
1
2
3
−5
4
0
1
ρ
2
3
4
ρ
Figura 4.2: Potenziale V (ρ) − E e dello spazio delle fasi per α = 0.5, ρb = 1 e
v = 0.25
in ρ = ρb valore di fondo per |θ| → ±∞, uno zero per ρ = ρm e uno zero per
ρ = ρa . La forma del potenziale cercata deve essere del tipo
Q(ρ) = 4P (ρ)(ρ − ρb )2 ,
(4.52)
con P (ρ) = aρ2 + bρ + c. Dal confronto con la (4.48) si ricavano i coefficienti
della forma quadratica
α
a=− ,
3
b=1−
2α
ρb
3
c=−
C02
= −v 2 ,
ρ20
(4.53)
quindi, le espressioni degli zeri ρm e ρa
ρm = −
−(3 − 2αρb ) +
(3 − 2αρb )2 − 12v 2 α
,
2α
#
(3 − 2αρb )2 − 12v 2 α
,
2α
e Q(ρ) può essere espresso in forma più compatta come
ρa =
(3 − 2αρb ) +
#
Q(ρ) = 4(ρ − ρa )(ρ − ρm )(ρ − ρ0 )2 .
(4.54)
(4.55)
(4.56)
L’espressione analitica della coppia di solitoni dark-antidark si ricava tramite
integrazione per quadratura della (4.48)
"
"
dρ
dθ =
(4.57)
Q(ρ)
in termini delle radici reali e1 < e2 < e3 < e4 del polinomio Q(ρ). La soluzione
della (4.57) è definita in due zone distinte: la prima per e1 < ρ < e2 e la seconda
per e3 < ρ < e4 . L’imposizione della condizione di periodicità ρ(θ) = ρ(θ+2π/k)
definisce il numero d’onda per le due soluzioni
* " e2
+−1
* " e4
+−1
dρ
dρ
#
#
k1 = π
,
k2 = π
(4.58)
Q(ρ)
Q(ρ)
e1
e3
45
4.2. SOLITONI DARK E ANTIDARK
Nella condizione di solitone i due zeri centrali collassano formando un doppio
zero del potenziale per e2 = e3 = ρb , e4 e e1 diventano rispettivamente ρm e ρa .
I profili di intensità ρ(θ) dei due solitoni si ottengono integrando la (4.57) tra
ρm e ρ, ottenendo il solitone dark, e integrando tra ρ e ρa che solitone antidark.
" ρ
" ρa
dρ
dρ
#
#
θ=
,
θ=
.
(4.59)
Q(ρ)
Q(ρ)
ρm
ρ
L’integrazione delle (4.59) è riportata in Appendice C, e la forma dei due solitoni
è rispettivamente
ρd =
dove
A=
ρm + rρa tanh2 Aθ
,
1 + r tanh2 Aθ
.
ρad =
α(ρa − ρb )(ρb − ρm )
,
3
r=
ρa + r" ρm tanh2 Aθ
,
1 + r" tanh2 Aθ
ρb − ρm
ρa − ρb
r" =
(4.60)
ρa − ρb
ρb − ρm
(4.61)
Sono riportate in figura 4.3 e 4.4 le forma del solitone dark e antidark per diversi
valori di α. Una volta note le soluzioni ρ(θ) si possono calcolare gli invarianti
1.2
1
α=0.2
α=0.5
d
ρ (θ)
0.8
α=0.8
0.6
0.4
0.2
0
−8
−6
−4
−2
0
θ
2
4
6
8
Figura 4.3: Profilo di intensità del solitone dark per ρb = 1 e v = 0.25
del moto
*
+
ρb
(u∗x u − ux u∗ ) 1 −
dx,
ρ
−∞
<
=
"
" ρ(x)
i +∞ |ux |2
H=
+
[g(ρb ) − g(ρ)] dρ dx,
2 −∞
2
ρb
M=
i
2
"
+∞
(4.62)
(4.63)
46
4.2. SOLITONI DARK E ANTIDARK
12
α=0.2
10
ρad(θ)
8
6
4
α=0.5
2
α=0.8
−8
−6
−4
−2
0
θ
2
4
6
8
Figura 4.4: Profilo di intensità del solitone antidark per ρb = 1 e v = 0.25
e la potenza, o il numero di particelle, Pa e Pd
" +∞
"
7 2
8
Pa =
|u| − ρ0 dx, ,
Pa =
−∞
+∞
−∞
7
8
ρ0 − |u|2 dx.
(4.64)
Le espressioni (4.60) permettono di ottenere, integrando la (4.42), la fase
non lineare associata alle due famiglie di solitoni. Si ottiene per il solitone dark
e antidark rispettivamente
*.
+
v
rρa
−
φd = − s tan 1
tanh [A(θ − θ0 )] + φ0
(4.65)
A
ρm
*.
+
v
rρm
−
φad = − s tan 1
tanh A(θ − θ0 )] + φ0
(4.66)
A
ρa
9
b −ρm )
dove s ≡ (ρa −ρρba)(ρ
. Dalle eq (4.65)-(4.66) si può calcolare il salto di fase
ρm
∆Φ = Φ(+∞) − Φ(−∞) lungo il solitone.
L’analisi dell’equilibrio del punto d’equilibrio ρ = ρb , ρ̇ = 0 mostra che esiste un
valore critico di α per cui le soluzioni solitoniche non sono più valide. Imponendo
∂Q(ρ)2 /∂ρ2 = 0 per ρ = ρb si ottiene αc :
αc =
ρb − v 2
ρ2b
(4.67)
É mostrato in figura 4.5 come varia lo spazio delle fasi quando α varia tra
zero e il valore critico αc . Prendendo l’intensità di fondo ρb = 1 e la velocità
47
4.2. SOLITONI DARK E ANTIDARK
2
a
0.1
c
V(ρ)−E
−2
e
0.05
0
0
V(ρ)−E
V(ρ)−E
2
−2
0
−0.05
−4
−4
0
2
4
0
2
ρ
4
4
b
−2
0
−2
2
ρ
4
−4
0.5
ρ
1
1.5
f
1
ρθ
0
0
0
2
d
2
ρθ
ρθ
2
−4
−0.1
4
ρ
0
−1
0
2
4
ρ
−2
0
0.5
ρ
1
Figura 4.5: Potenziale e spazio delle fasi per v = 0 e ρb = 1, e differenti valori di
α: (a-b) α = 0; (c-d) α = 0.5; (e-f) α = 0.9, dove la branca del solitone antidark
diventa estremamente piccola
v = 0, in accordo con la (4.54), si ha un solitone black (ρm = 0) per tutti i
valori di α. Il solitone antidark invece ha una intensità di picco ρa = 3/α − 2ρb
con una forte dipendenza da α che sarà spiegata in seguito. Nel limite α = 0 il
termine quintico nella (4.1) scompare e il modello decade nella NLS defocusing.
In questo limite il funzionale V (ρ) − E è cubico e dalle separatrici si dirama
una sola branca corrispondente alla soluzione di solitone dark della NLS (ρ =
tanh2 [x]). Invece, quando ρ > 1 il moto non è localizzato e non c’è coesistenza
di solitoni. Per 0 < α < αc si inverte il comportamento a θ = ∞ formando
la seconda ”buca” nel potenziale V (ρ) − E e il doppio loop della coppia darkantidark. Il solitone antidark è estremamente piccato per piccoli valori di α e
diventa progressivamente meno intenso quando α aumenta. Quando α → αc il
solitone antidark diventa estremamente debole fino a diventare un onda lineare
per α = αc con ρa → ρb .
Il legame dei solitoni dark-antidark con la velocità v si analizza rappresentando graficamente ρm e ρa in funzione di v per diversi valori di α (Fig. (4.6)).
Si può notare come il picco di intensità ρa resta costante per valori piccoli della
velocità v e decresce da infinito (α = 0) fino a ρa = ρb per α = αc quando il
solitone scompare. L’intensità di minimo ρm del solitone dark invece cresce con
la velocità formando un solitone diventa grey con ampiezza a = ρb −#
ρm , quindi,
la darkness del solitone decresce con v fino al valore minimo vc = ρb − αρ2b .
Da notare che il solitone della NLS (α = 0) è sempre meno intenso rispetto al
corrispondente della CQNLS con la stessa velocità.
1.5
48
4.2. SOLITONI DARK E ANTIDARK
1
0.9
γ=0.01
γ=0.6
0.8
γ=0.3
2
10
0.7
γ=0
γ=0.8
a
γ=0.1
0.5
ρ
ρm
0.6
0.4
1
10
γ=0.3
γ=0.9
0.3
0.2
γ=0.6
0.1
γ=0.9
0
0
0.1
0.2
0.3
0.4
0.5
0.6
velocità c
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
velocità c
0.7
0.8
0.9
1
Figura 4.6: (a) darkness del solitone dark e (b) massimo del solitone antidark
in funzione della velocità v per differenti valori di α e con ρb = 1.
4.2.1
Bistabilità e instabilità
La bistabilità nei sistemi solitonici implica l’esistenza di differenti soluzioni
aventi lo stesso valore di un invariante del moto (es., la potenza), al variare dei
parametri interni tipicamente ρb e v e la costante di propagazione non lineare
β = g(ρb ). Inoltre, Hermann [16] mostrò che i solitoni dark possono esibire
una bistabilità caratterizzata da differenti soluzioni aventi la stessa larghezza a
metà altezza (FWHM), sebbene in possesso di differenti ampiezze e, in generale,
differenti invarianti del moto. Questo secondo tipo di bistabilità fu introdotta
per distinguere dalla precedente definizione.
La coesistenza della coppia solitone dark-antidark costituisce un meccanismo
di bistabilità in quanto coesistono due differenti soluzioni con gli stessi parametri interni (ρb e v) e differenti invarianti M , H, P . Tuttavia, queste famiglie di
soluzioni possono essere bistabili anche secondo la definizione di Gantz e Herrmann: soluzioni con la stessa FWHM coesistono anche con differenti valori di
intensità di fondo ρb eα. L’espressione generale della FWHM, è
θF W HM =
#
2
tanh−1 f (ρ0 , ρm , ρa ),
w
(4.68)
ρb −ρm
a −ρb
dove f (ρ0 , ρm , ρa ) = 2ρaρ−ρ
per i solitoni dark e f (ρ0 , ρm , ρa ) = ρa +ρ
b −ρm
b −2ρm
per i solitoni antidark. I risultati ottenuti rappresentando θF W HM nel piano
(α, ρb ), in figura (4.7), mostrano che il solitone dark è bistabile per diversi
valori di α e tale range di valori di α si riduce notevolmente all’aumentare della
velocità. Al contrario, il solitone antidark non è mai bistabile secondo la seconda
definizione poiché non si osservano ripiegamenti delle curve di livello di θF W HM
per nessuna velocità v.
La stabilità della coppia di solitoni si può studiare applicando il criterio di
Barashenkov [17], ovvero analizzando le proprietà della derivata dell’invariante
del moto M (4.62) rispetto alla velocità v. La condizione marginale
∂M
= 0,
∂v
(4.69)
separa le soluzioni stabili (∂v M < 0) d quelle instabili (∂v M > 0). Il calcolo
di M per il solitone dark mostra che tale funzione è sempre negativa, quindi
l’intera famiglia di solitoni dark è stabile nel suo dominio di esistenza.
4.2. SOLITONI DARK E ANTIDARK
49
Figura 4.7: Curve di livello di FWHM nello spazio dei parametri α, ρb .(a)
solitone dark, (b) solitone antidark per v = 0.1. (c-d) solitone dark con v = 0.4
e v = 0.8. La linea tratteggiata rappresenta la condizione marginale αc .
Viceversa, i solitoni antidark mostrano un comportamento completamente
differente. Per piccoli valori di α il solitone antidark è sempre instabile su tutto
il dominio di esistenza. Per α sufficientemente grande il momento M (v) cambia
la sua pendenza vicino al cut-off vc invertendo la stabilità e formando solitoni
stabili.
La dinamica del solitone antidark, ad alti valori di α, per tre differenti velocità è riportata in figura 4.8. Nella 4.8(a) per v = 0.5, il solitone si propaga in
un fascio stabile (Fig 4.8b). Quando la velocità diminuisce, il solitone decade
in un differente solitone di maggiore velocità e stabile (Fig 4.8b). Infine, per
velocità piccole, il solitone decade in una coppia di solitoni antidark stabili (Fig
4.8c) che viaggiano con velocità opposte.
Al contrario, per piccoli valori di α non esiste soluzione stabile e il solitone
antidark decade in uno shock dispersivo (Fig 4.9). Per entrambe le direzioni, il
”trailing edge” è definito dal solitone. Il ”leading edge” è legato a onde piane che
aumentano la loro velocità e diventano sempre meno intense fino all’intensità di
fondo ρb = 1. Questo comportamento è analogo a quanto si può osservare nella
NLS (α = 0), in regime semiclassico, quando decade un impulso gaussiano su un
piedistallo [18]. Si sottolinea che in questo caso lo shock si forma direttamente
da un solitone instabile.
50
4.2. SOLITONI DARK E ANTIDARK
Figura 4.8: Dinamica del solitone antidark con ρb = 1, α = 0.7 per differenti
valori della velocità iniziale: (a) v=0.5, propagazione stabile; (b) v=0.4, decadimento in un solitone stabile; (c) decadimento in una coppia di solitoni antidark,
meno intensi, con velocità opposte
4.2.2
Parametrizzazione della soluzione di traslazione negli shock dispersivi
La soluzione del modello (4.60) dipende dai parametri interni v e ρb e un
parametro esterno α. La scelta di tali parametri è completamente arbitraria
e dipende esclusivamente da cosa si vuole studiare nel modello. Negli shock
dispersivi si preferisce usare ρm , l’intensità di minimo del solitone dark e ricavare
l’espressione della velocità v in funzione di essa. Questa parametrizzazione si
può ottiene imponendo una diversa espressione analitica del potenziale Q(ρ)
nel limite solitonico. Si consideri dalla soluzione del sistema in un funzione
dell’intensità ρ
* +2
dρ
= Q(ρ),
(4.70)
dθ
4
3
2
con Q(ρ) = − 4α
3 ρ + 4ρ + a2 ρ + a3 ρ + a4 , dove a2 , a3 , a4 sono
a2 = −8k0 ,
a3 = 8E0 ,
a4 = −4C02 .
(4.71)
La condizione di solitone porta al fenomeno collasso degli zeri formando un
doppio zero della funzione Q(ρ) in ρb e altri due zeri per ρ = ρm e ρ = ρa . Dato
che ρm in questa parametrizzazione è un parametro interno del sistema si cerca
un potenziale Q(ρ) del tipo
Q(ρ) = 4(ρ − ρb )2 (ρ − ρm )(Aρ − B),
(4.72)
4.2. SOLITONI DARK E ANTIDARK
51
Figura 4.9: Decadimento di un solitone antidark instabile (ρb = 1, v = 0, α =
0.1) in due shock dispersivi simmetrici: (a) grafico di intensità; (b) istantanea
dell’ intensità e della fase per z = 5.
Sviluppando la (4.72) e confrontando con la (4.70) si determina la forma dei
parametri A e B
,
α
α
A=− ,
B = 1 − (2ρb + ρm ) ,
(4.73)
3
3
quindi, il terzo zero ρa = −B/A, intensità di picco solitone antidark è
ρa =
3
− (2ρb + ρm ).
α
(4.74)
Confrontando la (4.72) con la (4.56) si determina la relazione tra v e ρm
/
0
7
a 8
v 2 = ρm 1 − αρb 1 −
.
(4.75)
3ρb
dove a = ρb − ρm . Il motivo per cui si preferisce questa parametrizzazione è
legato al cambiamento di variabile delle equazioni di modulazione, menzionato
nel capitolo 3, nella formulazione non integrabile della costruzione di GurevichPitaevskii. L’ampiezza a del solitone dark è, con questa parametrizzazione,
una espressione semplice dei soli parametri interni ρb e ρm e non contiene il
parametro esterno α. Si introduce, inoltre, un ulteriore parametro che sarà
52
4.3. SHOCK DISPERSIVO
utilizzato nello studio dello shock dispersivo, κ. Quando |θ| →∞ . si ha che
l’intensità del solitone decade come
ρb − ρ ∝ exp(−κ|θ|)
(4.76)
Si può trovare la relazione tra κ e gli altri parametri analizzando la (4.48) nel
limite |θ| → ∞ dove ρ → ρb . Espandendo in serie il potenziale Q(ρ) per piccoli
valori di ρ" = ρb − ρ si ha
6
6
dQ(ρ) 66 " 1 d2 Q(ρ) 66
Q(ρ) = Q(ρb ) +
ρ +
(ρ" )2 ,
(4.77)
dρ 6ρb
2 dρ2 6ρb
dove i primi due termini sono nulli per le proprietà legate al solitone. Mettendo
in relazione la (4.77) con la (4.48) e la (4.76) si ottiene
6
1 d2 Q(ρ) 66
2 " 2
κ (ρ ) =
(ρ" )2
(4.78)
2 dρ2 6ρb
quindi,
κ=
4.3
4
6
#
1 d2 Q(ρ) 66
= 2 f (ρb ) − v 2
2 dρ2 6ρb
(4.79)
Shock dispersivo
La soluzione periodica della CQNLS (4.60) espressa in termini della fase
veloce θ è caratterizzata da tre costanti della riduzione (4.71) e la velocità v
(4.75). In una teoria di modulazione questi parametri diventano lentamente
variabili rispetto ad x e z. Nella teoria di Whitham le equazioni di modulazione
si ricavano dalle leggi di conservazioni dell’equazione dispersiva mediata su una
oscillazione veloce rispetto alla fase θ:
∂
∂
(ρ) +
(ρu) = 0,
∂z
∂x
/ 2
0
∂
∂ u
ρxx
(ρx )2
(u) +
+ f (ρ) +
−
= 0,
∂z
∂x 2
4ρ
8ρ2
/
0
∂
∂
ρ 1
(ρu) +
ρu2 + f (ρ) − (ρ log(ρ))xx = 0,
∂z
∂x
2 2
(4.80a)
(4.80b)
(4.80c)
che unita con la legge di conservazione del numero d’onda k compongono le
equazioni di modulazione. Il sistema di Whitham per la GNLS può essere
espresso nella forma:
:
;
:
;
P i (a2 , a3 , a4 , v) z + Qi (a2 , a3 , a4 , v) x = 0,
i = 1, 2, 3,
(4.81)
k(a2 , a3 , a4 , v)z + ω(a2 , a3 , a4 , v)x = 0.
(4.82)
Dove P1 = ρ, P2 = u e P3 = ρu sono le densità conservate della GNLS e
Qi i corrispondenti flussi. Le medie sono calcolate sulla famiglia di soluzioni
periodiche (4.58)
"
k e2 f (a2 , a3 , a4 , v; ρ)
#
f (a2 , a3 , a4 , v) =
dρ
(4.83)
π e1
Q(ρ)
53
4.3. SHOCK DISPERSIVO
La descrizione analitica dello shock dispersivo prevede di risolvere il sistema
(4.81), (4.82), con condizioni al contorno tali da garantire continuità del sistema
con il limite non dispersivo ai bordi x± (z) dello shock dispersivo
x = x+ (z) :
a = 0,
ρ = ρ+ ,
u = u+ ,
(4.84)
x = x− (z) :
k = 0,
ρ = ρ− ,
u = u− ,
(4.85)
La dipendenza di ρ, u, k, a rispetto a a2 , a3 , a4 e v. è definita nella (4.83) e
dalle formule nella sezione $4.3. Le coppie (ρ− , u− ) e (ρ+ , u+ ) rappresentano
le soluzioni del limite non dispersivo (4.6) valutati rispettivamente al ”trailing
edge” e al ”leading edge”’. I bordi x± (t) dello shock dispersivo rappresentano
le condizioni al contorno libere che contengono le condizioni cinematiche:
dx+
= v0 (ρ+ , u+ , k + ) ,
dz
dx−
= vs (ρ− , u− , a− ) ,
dz
(4.86)
dove v0 (ρ+ , u+ , k) = ∂ω0 /∂k è la velocità di gruppo di un pacchetto d’onda
lineare con numero d’onda k che si propaga su intensità di fondo ρ+ e fase u+
e vs (ρ− , u− , a) è la velocità del solitone con ampiezza a e variabili di fondo ρ−
e u− . Ovviamente, i valori di k + al ”leading edge”, l’ampiezza a− al ”trailing
edge” e le velocità dei bordi x± (t) devono essere determinati come previsto dal
problema di Gurevich-Pitaevskii.
Si determineranno ora le velocità s± dei bordi, secondo la riformulazione del
problema di Gurevich-Pitaevskii per modelli non integrabili esposta nel Cap. 3,
utilizzando le semplificazioni che le equazioni di modulazione subiscono quando
a = 0 o k = 0. La costruzione verrà applicata esclusivamente al solitone dark.
Al ”leading edge” il calcolo della velocità s+ è lo stesso sia per il solitone dark
che per il solitone antidark. Al ”trailing edge” invece, la soluzione per il solitone
antidark, secondo gli argomenti esposti nella sezione 4.2.3, si ottiene cambiando
parametro interno ρm → ρa e mandando a → −a.
4.3.1
Leading Edge
Al ”leading edge”, quando x = x+ (z), l’ampiezza di oscillazione a → 0.
Come per la KdV, si deve adottare la procedura di riduzione della teoria di
Whitham esposta nel Cap. 3. Se a = 0 la media di una osservabile rispetto a
due variabili è pari all’osservabile calcolata per le due variabili mediate
F (ρ, u) = F (ρ, u).
(4.87)
La riduzione delle equazioni di Whitham deve essere consistente con le condizioni
al contorno dello shock dispersivo quando a = 0. Al ”leading edge” le soluzioni devono raccordarsi con il limite dispersivo e le equazioni di modulazione si
riducono a
ρz + (ρ u)x = 0,
(4.88)
uz + u ux + f " (ρ)ρx = 0.
Le medesime equazioni possono essere ottenute direttamente dal sistema di Whitham imponendo a = e2 − e1 → 0 nelle (4.81). Inoltre, la medesima riduzione
deve riflettersi anche nella legge di conservazione del numero d’onda, che generalmente in un modello integrabile permette la chiusura del sistema. La chiave della
riduzione sta nella scelta del cambiamento di variabili menzionato nel paragrafo
54
4.3. SHOCK DISPERSIVO
$3.1. Per un sistema del tipo GNLS il cambiamento è (a2 , a3 , a4 , v) → (ρ, u, k, a).
La legge di conservazione del numero d’onda espressa nel nuovo set di variabili
è
kz + [ω(ρ, u, k, a)]x = 0,
con
ω = vk,
(4.89)
che al ”leading edge” diventa:
dove
kz + [ω0 (ρ, u, k)]x = 0,
(4.90)
.
/
0
k2
2
ω0 (ρ, u, k) = k u +
+ (ρ − αρ ) ,
4
(4.91)
u = R(ρ) − R(1)
(4.92)
è la relazione di dispersione di un pacchetto d’onda lineare. Le equazioni (4.88)
e (4.90) compongono un sistema chiuso che rappresenta la riduzione esatta per
ampiezza nulla delle equazioni di Whitham.
Nel sistema ridotto si riconosce l’equivalente del sistema idrodinamico con l’aggiunta della conservazione del numero d’onda. I due sistemi sono disaccoppiati
quindi possono essere risolti indipendentemente per ρ(x, z) e u(x, z). Dato che
i valori di ρ e u sono soggetti alle condizioni al contorno per a = 0, si può ammettere che la diagonalizzazione del sistema idrodinamico è eseguibile tramite
la ”simple-wave”. Imponendo la relazione su u
e sostituendo la (4.92) nella seconda del sistema (4.88) si ha
.
.
*
+
1 − αρ
1 − αρ
ρt +
R(ρ) − R(1) ρx + (1 − αρ)ρx = 0,
ρ
ρ
(4.93)
quindi, l’espressione in forma diagonale della equazione di modulazione
ρz + V+ (ρ)ρx = 0,
(4.94)
9
(4.95)
dove
V+ (ρ) = R(ρ) − R(1) +
ρ − αρ2 .
Invece, usando la definizione di u nella (4.90) si ottiene una seconda equazione
differenziale dipendente dalle sole variabili (ρ, u)
kz + [Ω(ρ, k)]x = 0,
(4.96)
con
.
/
0
k2
Ω(ρ, k) = ω0 (ρ, u(ρ), k) = k R(ρ) − R(1) +
+ ρ − αρ2 .
4
(4.97)
Riconosciamo nelle equazioni (4.94) e (4.97) la stessa riduzione ottenuta per la
KdV (3.24) con le velocità definite come
V+ (ρ),
∂Ω(ρ, k)
∂k
(4.98)
55
4.3. SHOCK DISPERSIVO
che definiscono due famiglie di curve caratteristiche:
dx
= V+ (ρ)
dz
(4.99)
e
dx
∂Ω(ρ, k)
=
(4.100)
dz
∂k
La famiglia definita dalla (4.99) dipende unicamente da ρ ed è completamente
determinata dalla simple wave del limite dispersivo della GNLS, l’equivalente
dell’equazione di Hopf nella KdV. Questa famiglia di soluzioni collega la zona interna in cui termina lo shock dispersivo con la zona in cui vale il limite
idrodinamico. La seconda famiglia di soluzioni, invece, dipende da ρ e k. La
definizione delle velocità s+ del ”leading edge” x = x+ (z) in accordo con la
(4.97) mostra che il ”leading edge” dello shock dispersivo è rappresentato dalla
caratteristica (4.100).
Come per la KdV, sostituendo k = k(ρ) nelle equazioni (4.94) e (4.96) si ottiene
la relazione che lega le due velocità caratteristiche
dk
∂Ω/∂ρ
=
dρ
V+ − ∂Ω/∂k
(4.101)
L’equazione differenziale in k si risolve con la condizione iniziale k(ρ− ) = 0 che
equivale alla formazione del solitone al ”traling edge”. É utile definire una nuova
variabile γ
4
k2
γ = 1+
(4.102)
4(ρ − αρ2 )
Il numeratore e il denominatore della (4.101) espresse in funzione di γ diventano
/
0
∂Ω
k
#
=
γ + (1 − 2αρ)(γ + 1))
(4.103)
∂ρ
2γ ρ − αρ2
e
∂Ω
V+ −
=
∂k
#
/
0
ρ − αρ2
(1 − γ)(2γ + 1) .
γ
(4.104)
Il differenziale dk/dρ espresso in funzione di γ si ricava invertendo la (4.102)
9
k 2 = 4(γ 2 − 1)( ρ − αρ2 ),
(4.105)
e derivando rispetto a ρ
2k
dk
dγ
= 4(γ 2 − 1)(1 − 2αρ) + 8γ(ρ − αρ2 )
dρ
dρ
(4.106)
Sostituendo le (4.103),(4.104) e(4.106) nella (4.101) si ottiene una nuova equazione differenziale in γ
dγ
(1 + γ)[1 + 2γ − 4γαρ]
=−
dρ
2(2γ + 1)(ρ − αρ2 )
(4.107)
con la condizione iniziale γ(ρ− ) = 1, dove ρ− è determinata in funzione della
discontinuità iniziale ρ0 . Se è nota la soluzione della (4.107), quindi, l’espressione
56
4.3. SHOCK DISPERSIVO
della γ(ρ), si risale al numero d’onda k + al ”leading edge” , dove ρ = 1, tramite
la (4.105)
#
k + = k(1) = 2(1 − α) γ 2 (1) − 1.
(4.108)
La velocità di propagazione s+ del ”leading edge” è
*
+
dx+
∂Ω(1, k + )
1
+
s =
=
= (1 − α) 2γ(1) −
dz
∂k
γ(1)
(4.109)
L’equazione (4.107) non è risolvibile analiticamente, tuttavia è possibile trovare
delle correzioni a partire dal caso limite α = 0 i cui la la GNLS si riduce alla
NLS. Infatti, se α = 0, la (4.107) diventa
dγ
1+γ
=−
,
dρ
2ρ
la cui soluzione è
con γ(ρ− ) = 1,
4
γ(ρ) = 2
ρ−
− 1,
ρ
(4.110)
(4.111)
e la velocità s+ diventa
s+ =
#
8ρ− − 8 ρ− + 1
#
2 ρ− − 1
(4.112)
in completo accordo con i risultati noti per la NLS [10]. Per valori di α $ 1, si
possono trovare delle correzioni alla (4.111). Introducendo
γ(ρ) = γ0 + γ1 ,
(4.113)
dove γ0 soddisfa la (4.111) e γ1 è una correzione di ordine α. Espandendo in
serie la (4.107), trascurando i termini di ordine α2 , si arriva all’equazione per
γ1
4
#
4 ρ− /ρ − 3 ρ−
dγ1
γ1
=− + #
α,
(4.114)
dρ
2ρ 4 ρ− /ρ − 1
ρ
che integrata con la condizione iniziale γ1 (ρ− ) = 0 da come soluzione
!
+
#
#
# / * 16
# 0>
1
−1
−
−
−
−
γ1 (1) = 2 ρ α 3+16ρ −19 ρ +32 ρ ln
− #
+tanh 4 ρ− .
9
9 ρ−
(4.115)
Sostituendo l’espressione di γ1 (1) nella (4.109) si ottiene l’espressione esplicita
della velocità s+
#
/
0
8ρ− − 8 ρ− + 1
1
+
#
#
s ≈
(1 − α) + 2 −
γ1 (1),
(4.116)
2 ρ− − 1
(2 ρ− − 1)2
a meno di termine dell’ordine α2 .
57
4.3. SHOCK DISPERSIVO
4.3.2
Trailing Edge
Quando x = x− (z), al ”trailing edge”, si genera un treno di solitoni dark
debolmente interagenti che si propagano su un fondo ρ, u lentamente variabile.
Dato che k → 0 per x → x− , si deve analizzare il limite solitonico del sistema di
modulazione (4.81). In questo limite, la lunghezza d’onda L tende a infinito, il
contributo dei solitoni nei flussi mediati ρ e u tendono a zero, quindi, come per
il ”leading edge” si ha che F (ρ, u) = F (ρ, u). La differenza fondamentale è che
il limite k → 0 è singolare e richiede un approccio diverso. Si vuole utilizzare
il metodo dell’onda coniugata esposto nel precedente capitolo. Ovvero si vuole
risolvere il problema solitonico del sistema originale tramite l’analisi del limite
armonico del sistema coniugato.
Il primo passo è introdurre la trasformazione di coordinate
x → ix̃
z → it̃
(4.117)
θ → iθ̃.
nella (4.4). Si ottiene un sistema identico all’originale a meno di un segno nei
termine dispersivi
ρt + (ρu)x = 0
$
xx
+
ut + uux + f " (ρ)ρx = − ρ4ρ
ρ2x
8ρ2
%
(4.118)
x
Imponendo le trasformazioni (4.117) alla (4.70) si ottiene la soluzione associata
al sistema coniugato
* +2
dρ
= −Q(ρ) ,
(4.119)
dθ?
Il segno meno determina un ribaltamento del potenziale Q(ρ) il cui effetto è invertire la zona solitonica con quella armonica. Dando le condizioni di periodicità
ρ(θ̃) = ρ(θ̃ + 2π/k̃) si definisce il numero d’onda coniugato
?
k=π
@"
e3
e2
#
dρ
−Q(ρ)
A−1
.
(4.120)
Nel limite solitonico, e2 → e3 , si può espandere Q(ρ) intorno al punto ρ = e2 =
e3 e la (4.119) diventa la legge di conservazione dell’energia per l’oscillatore
armonico:
6
* +2
1 dρ
1 d2 Q 66
+
(ρ − ρ)2 = Q(ρ).
2 dθ?
4 dρ2 6ρ
e dal confronto con la (4.79) si mostra che in questo limite
4
6
1 d2 Q 66
?
k=
= κ,
2 dρ2 6ρ
(4.121)
La legge di dispersione dell’onda coniugata si ottiene linearizzando le equazioni
(4.118) attorno a valori di fondo ρ̄ e ū:
.
*
+
κ2
ω
? (ρ, u, k̃) = ω
? (ρ, u, κ) = κ u + f (ρ̄) −
,
(4.122)
4
58
4.3. SHOCK DISPERSIVO
che, dopo la sostituzione di u nel limite non dispersivo (4.88) e f (ρ̄), assume la
forma
.
/
0
κ2
2
?
Ω(ρ, κ) = ω
?0 (ρ, u(ρ), k) = κ R(ρ) − R(1) + ρ − αρ −
.
(4.123)
4
La trasformazione (4.117) mappa un’onda armonica exp[i(kx − ωz)] in un’onda
solitonica exp[±κ(x − vs z)], con una legge di dispersione solitonica ω
? = vs κ. Nel
limite solitonico dove k̃ = κ, la le legge di dispersione lineare ω (4.91) e la legge
di dispersione solitonica ω
? (4.122) sono legate dalla trasformazione
i?
ω (κ) = ω(iκ).
(4.124)
Ora che è stata dimostrata l’equivalenza dei due sistemi in questo limite si
può procedere allo stesso modo del ”leading edge” con una ridefinizione delle
variabili. Tramite considerazioni identiche a quanto esposto per la KdV, si può
sostituire l’ampiezza a con il numero d’onda coniugato κ e il numero d’onda
originale con
Λ ≡ k/k̃
(4.125)
Il sistema di Whitham nel limite non dispersivo dell’onda coniugata è identico
a quello originale (4.88) e si ha nuovamente che le equazioni ridotte devono
raccordarsi con la simple-wave
ρz + V+ (ρ)ρx = 0,
(4.126)
con V+ definita nella (4.95). Il passaggio successivo è eseguire il limite k → 0
nella legge di conservazione (4.82). Introducendo il cambiamento di coordinate
(a2 , a3 , a4 , v) → (ρ̄, ū, k̃, Λ) e operando la sostituzione k = Λκ si ha
k̃Λz + ω̃Λx + Λ(k̃z + ω̃x ) = 0 ,
(4.127)
con ω
? = v k̃. Si ottiene un’equazione identica alla (3.46) a meno di una ridefinizione di ω. L’espansione per Λ $ 1 conduce all’espressione dell’equazione
caratteristica per i termini dominanti.
per cui
? ∂Λ
∂Λ Ω
+
= 0,
∂z
κ ∂x
(4.128)
dx
Ω̃(ρ̄, κ)
=
,
(4.129)
dz
κ
dove Λ0 $ 1 è una costante. In particolare, quando Λ0 = 0 la caratteristica (4.129) definisce completamente il ”traling edge”. Considerando la (4.127)
?
? si
calcolata lungo la curve caratteristiche dx/dz = Ω/κ
dove ?
k = κ, ω
? = Ω
ricava
κz + Ω x = 0
(4.130)
Λ =Λ 0
per
Dato che ρ e κ non possono essere specificate indipendentemente lungo una
caratteristica deve esistere una relazione locale κ(ρ) consistente con il sistema
(4.126), (4.130). Sostituendo κ = κ(ρ) nella (4.130) e usando la (4.126) si ottiene
?
dκ
∂ Ω/∂ρ
=
.
?
dρ
V+ − ∂ Ω/∂κ
(4.131)
59
4.3. SHOCK DISPERSIVO
Le velocità caratteristiche al ”trailing edge” sono
V+ ,
Ω̃
.
κ
(4.132)
dove il primo termine V+ rappresenta la velocità del limite non dispersivo, la
simple wave, e il secondo termine è la velocità che connette la soluzione interna
con quella esterna. La condizione iniziale per l’equazione differenziale ordinaria (4.132) deve essere applicabile anche nel caso di intensità nulla dello shock
dispersivo che è una condizione globale (4.85), corrispondente ai valori iniziali
di ρ− = ρ+ = 1. Si richiede quindi κ(1) = 0. In accordo con le condizioni
cinematiche (4.86), la velocità del bordo è pari alla velocità del solitone, e si ha
s− =
? − , κ− )
Ω(ρ
dx−
=
,
dz
κ−
dove κ− = κ(ρ− ). Introducendo la nuova variabile
4
κ2
γ
?= 1−
4(ρ − αρ2 )
(4.133)
(4.134)
nella (4.131) e tramite operazioni analoghe effettuate al ”leading edge” si ottiene
l’equazione differenziale ordinaria
dγ̃
(1 + γ̃)[1 + 2γ̃ − 4γ̃αρ]
=−
dρ
2(2γ̃ + 1)(ρ − αρ2 )
con la condizione iniziale
γ
?(1) = 1.
(4.135)
(4.136)
Nota la funzione γ
?(ρ), la velocità al ”trailing edge” è determinata dalla (4.133)
unita con la (4.123)
#
s− = R(ρ− ) − R(1) + ρ − αρ2 γ
?.
(4.137)
L’ampiezza a = ρ− − ρm del solitone al ”trailing edge” in funzione del salto di
intensità ρ− e di γ
? allo shock dispersivo si ricava dalla (4.75) con v = s− − u− ,
−
ρb = ρ
*
+
/
*
+0
7 82
a
ρ− − α ρ−
γ
?2 = (ρ− − a) 1 − αρ− 1 − −
(4.138)
3ρ
Come per il ”leading edge”, quando α = 0 il sistema decade nella NLS di cui si
può scrivere la soluzione esplicitamente. Dalla (4.135) si ottiene
d?
γ
1+γ
?
=−
,
dρ
2ρ
(4.139)
la cui soluzione con le condizioni al contorno (4.136) è
2
γ
?(ρ) = √ − 1.
ρ
(4.140)
60
4.3. SHOCK DISPERSIVO
Mentre, le equazioni (4.137) e (4.138), unite con la (4.140), danno rispettivamente
#
s− = ρ−
(4.141)
e
#
a = 4( ρ− − 1)
(4.142)
in completo accordo con risultati noti [10]. Ancora, per piccoli valori di α si
possono trovare correzioni della (4.139) in forma esplicita. Introducendo γ
? =
γ
?0 + γ
?1 , dove γ
?0 è data dalla (4.140), γ
?1 deve soddisfare
√
d?
γ1
γ
?1
4/ ρ − 3 α
√ ,
=− + √
(4.143)
dρ
2ρ 4/ ρ − 1 ρ
che integrata con la condizione iniziale γ
?1 (1) = 0,
=C
B
<
#
#
−
−
−−1
4
−
ρ
ρ
α
ρ
−
1
γ
?1 (ρ− ) = #
+
+
,
ρ− − 1 + 64 ln
3
4
32
ρ−
(4.144)
e quindi
/
0
#
#
2 −# −
−
−
−
−
−
s ≈ ρ (1 + γ
?1 (ρ )) − (ρ
ρ − 1) + ρ (2 − ρ ) α.
(4.145)
3
Capitolo 5
Conclusioni
In questa tesi è stata esposta una trattazione completa per la risoluzione dello
shock dispersivo includendo un esempio pratico nel caso della Korteweg-de Vries
e dell’equazione non lineare di Schrödinger a nonlinearità quintica (CQNLS). Il
metodo di risoluzione esposto è la teoria di Whitham nelle sua formulazione per
sistemi integrabili e non integrabili.
A livello generale, la teoria di Whitham caratterizza completamente lo shock
dispersivo di un modello. In particolare si calcolano le velocità di propagazione
dei bordi dello shock, tramite le equazioni di modulazione dei parametri lentamente variabili dell’onda cnoidale, quest’ultima soluzione dell’equazione iniziale.
La costruzione di Gurevich-Pitaevskii richiede che le equazioni di modulazione,
o di Whitham, siano espresse in forma diagonale e quindi necessita del calcolo
degli invarianti di Riemann. Tuttavia, il calcolo degli invarianti rappresenta il
principale ostacolo della teoria perchè, escludendo casi particolari (KdV e NLS),
è di difficile soluzione.
La formulazione di Whitham per i sistemi non integrabili utilizza delle proprietà di raccordo delle leggi di conservazione con il limite non dispersivo del
modello riuscendo a ridurre il numero delle equazioni di modulazione. Con
questo approccio, la risoluzione dello shock dispersivo passa per il calcolo degli
invarianti di Riemann del sistema ridotto e non del sistema di Whitham generale. Si sottolinea che l’esistenza degli invarianti di Riemann del sistema ridotto
non influenza l’esistenza degli invarianti del sistema generale. Pur non calcolando la forma diagonale leggi di modulazione dell’onda cnoidale, tramite le
equazioni del sistema ridotto, diagonalizzate, si riesce a caratterizzare lo shock
dispersivo calcolandone i parametri principali.
In questa tesi è stata esposta, in maniera completa, l’applicazione della formulazione non integrabile per la CQNLS. Tale equazione ammette fenomeni di
shock, in regime semiclassico, legati dal decadimento di un solitone dark. Tale
esempio mostra chiaramente che nei sistemi bidirezionali si ottiene una riduzione
delle equazioni di Whitham a sistemi 2x2, sempre diagonalizzabili (Appendice
A). Questa riduzione ha permesso di calcolare i parametri principali dello shock,
in particolare, le velocità s± dei bordi e l’ampiezza del solitone al trailing edge
in funzione del salto di intensità iniziale.
61
62
Inoltre, per la CQNLS, si è mostrata l’esistenza di due differenti soluzioni,
la coppia di solitoni dark-antidark, legate da un meccanismo di bistabilità del
primo e del secondo tipo. In aggiunta è stato studiato un nuovo meccanismo di
instabilità che coinvolge la formazione di un DSW a seguito del decadimento di
un solitone antidark instabile.
Appendice A
Invarianti di Riemann
Gli invarianti di Riemann permettono di riscrivere in forma più agevole un
sistema di equazioni differenziali alle derivate parziali. Questo metodo è una
generalizzazione del metodo delle caratteristiche, ovvero, la ricerca di funzioni
che si mantengano costanti lungo le curve. In questo paragrafo si espone un
metodo per la determinazione degli invarianti in un sistema 2x2.
Si consideri il generico sistema quasi lineare con due variabili ρ(x, t) e u(x, t):
a1 ρt + a2 ρx + a3 ut + a4 ux = 0,
b1 ρt + b2 ρx + b3 ut + b4 ux = 0.
(A.1)
Si consideri la curva α nel piano (t, x) definita dalle espressioni
t = T (s),
x = X(s),
(A.2)
dove s è un parametro. Si definiscono le componenti del vettore tangente alla
curva α
dx
dt
= X " (s) = cos φ
= T " (s) = sin φ,
(A.3)
ds
ds
dove φ rappresenta l’angolo formato tra il vettore tangente e l’asse x. Supponendo che i valore iniziali ρ(x, t) e u(x, t) siano dati sulla curva α, le derivate
direzionali di u e ρ sono
ρ" ≡ dρ/ds = ρx cos φ + ρt sin φ,
u" ≡ du/ds = ux cos φ + ut sin φ.
(A.4)
Anche il sistema (A.1) deve essere soddisfatto, quindi ci sono quattro equazioni
per la determinazione delle derivate ρx , ρt , ut e ux . Dato che il dato iniziale è
dato sulla curva caratteristica, il sistema non ha soluzioni, quindi il determinante
deve essere nullo. Si ha che per date soluzioni ρ(x, t) e u(x, t) l’angolo φ deve
soddisfare la condizione
6
6 A1
6
6 B1
quindi,
6 6
A2 66 66 a1 cos φ − a2 sin φ
=
B2 6 6 b1 cos φ − b2 sin φ
6
a3 cos φ − a4 sin φ 66
= 0,
b3 cos φ − b4 sin φ 6
∆13 cos2 φ − (∆14 − ∆23 ) sin φ cos φ + ∆24 sin2 φ = 0,
63
(A.5)
(A.6)
64
dove
∆ij = ai bj − bi aj .
(A.7)
Se la (A.6) ha due soluzioni reali in termini di cot φ si possono considerare due
famiglie di curve caratteristiche corrispondenti alle soluzioni ρ(x, t) e u(x, t). In
questo caso il sistema (A.1) è detto iperbolico. Se si ha una sola soluzione il
sistema è parabolico, invece, per due soluzioni complesse il sistema è ellittico.
Supponendo che tali radici siano reali, vogliamo ricavare come le funzioni ρ e u
variano lungo le curve caratteristiche.
L’idea generale è di esprimere le (A.1), in termini delle derivate direzionali di ρ e
u lungo la curva. Dal sistema (A.1), sottraendo la prima equazione moltiplicata
per B1 con la seconda moltiplicata per A1 , si ottiene
−∆12
dρ
+ ux cos φ(∆23 − ∆13 cot φ) + ut sin φ(−∆14 + ∆24 tan φ) = 0. (A.8)
ds
Alcune manipolazioni della (A.6) unite con l’antisimmetria della (A.7) portano alle identità trigonometriche
−∆14 − ∆13 cot φ = −∆23 − ∆42 tan φ ≡ −K,
−∆32 − ∆13 cot φ = −∆14 − ∆42 tan φ ≡ −L.
(A.9)
Imponendo la seconda delle (A.9) nella (A.8) si ricava
∆12
dρ
du
+K
= 0.
ds
ds
(A.10)
Allo stesso modo, sottraendo la prima delle (A.2) moltiplicata per B2 con la
seconda moltiplicata per A2 e si ottiene
L
dρ
du
+ ∆34
= 0.
ds
ds
(A.11)
Se nessuna delle equazioni (A.10) e (A.11) ammettono la soluzione identica, esse
sono equivalenti. Inoltre, K e L dipendono da cot φ, quindi la scelta della radice
della (A.6) definisce la variazione lungo una delle curve caratteristiche C1 e C2 .
Prendendo, per esempio, la variazione di ρ e u definite nella (A.10) si ha
∆12
dρ
du
+ K(cot φ1 )
= 0,
ds1
ds1
∆12
dρ
du
+ K(cot φ2 )
= 0.
ds2
ds2
(A.12)
con s1 e s2 i parametri relativi rispettivamente alla famiglia C1 e C2 . Si aggiunge, inoltre, l’ipotesi che i coefficienti ai e bi non dipendano esplicitamente dalle
coordinate x e t, ma solo attraverso ρ(x, t) e u(x, t).
Se si riescono a trovare delle funzioni λ1,2 = λ(ρ, u) tali da soddisfare
γ∆12 = ∂λ1 /∂ρ,
β∆12 = ∂λ2 /∂ρ,
γK(cot φ1 ) = ∂λ1 /∂u,
βK(cot φ2 ) = ∂λ2 /∂u.
(A.13)
∂(β∆12 )
∂(βK(cot φ2 ))
=
∂u
∂ρ
(A.14)
e le condizioni di compatibilità
∂(γ∆12 )
∂(γK(cot φ1 ))
=
,
∂u
∂ρ
65
è possibile riscrivere la (A.10) nella forma
dλ1
= 0,
ds1
dλ2
= 0.
ds2
(A.15)
Le due variabili λ1,2 si mantengono costanti lungo le caratteristiche e sono chiamate invarianti di Riemann. Le due curve definiscono in accordo con la (A.3) i
due vettori tangenti (cos φ1,2 , sin φ1,2 ) nel piano (x, t) e le derivate (A.10) hanno
il significato di variazioni lungo le curve integrali definite dai campi vettoriali.
Quindi possono essere riscritte
∂λ1
∂λ1
cos φ1 +
sin φ1 = 0,
∂x
∂t
oppure
dove v1,2 = cot φ1,2
∂λ2
∂λ1
cos φ2 +
sin φ2 = 0,
∂x
∂t
∂λ1,2
∂λ1,2
+ v1,2
= 0,
∂t
∂x
sono le velocità caratteristiche.
(A.16)
(A.17)
Invarianti nella NLS quintica Partendo dal sistema idrodinamico derivante
dalla NLS quintica
ρt + (ρu)x = 0
(A.18)
ut + uux + f " (ρ)ρx = 0
dove f (ρ) = ρ + 2δ ρ2 e δ = ηγ, con η = ±1. Si ricavano i coefficienti ai , bj :
a1 = 1,
b1 = 0,
a2 = u,
b2 = f " (ρ),
a3 = 0,
b3 = 1,
a4 = ρ,
a4 = u.
(A.19)
I determinanti ∆ij = ai bj − bi aj sono:
∆13 = 1,
∆14 = ∆23 = u,
∆24 = u2 − ρf " (ρ).
(A.20)
Date le curve caratteristiche ρ(x, t) e u(x, t) soluzioni del problema (A.18), la
relazione tra l’angolo φ della tangente alla curva caratteristica e l’asse x è definita
dall’equazione (A.6), con i coefficienti espressi nella (A.19). Si ha quindi
cos2 φ−2u sin φ cos φ+(u2 −ρf " (ρ)) sin2 φ ≡ (cos φ−u sin φ)2 −ρf " (ρ) sin2 φ = 0.
(A.21)
Risolvendo la (A.21) in funzione di cot φ si ottengono le velocità caratteristiche
#
v± = cot φ± = u ± ρf " (ρ)
(A.22)
nel modello descritto dalla (A.18) le velocità caratteristiche assumono la forma:
#
V± = u ± ρ + δρ2 .
(A.23)
#
Il coefficiente K nella (A.12) è pari a K(cot φ± ) = ± ρ + δρ2 e ∆12 = 1 − δρ.
Le condizioni di compatibilità (A.14) diventano:
%
∂ $
∂ #
α(1 + δρ) =
(α ρ + δρ2 )
∂u
∂ρ
(A.24)
66
%
∂ $
∂ #
β(1 + δρ) = − (β ρ + δρ2 )
∂u
∂ρ
Si osserva che le equazioni (A.24) e (A.25) sono soddisfatte se
4
1 + δρ
γ=β=
.
ρ
(A.25)
(A.26)
Si ottiene dalle relazioni (A.13)
∂λ1
(1 + δρ)
=#
,
∂ρ
ρ + δρ2
e
∂λ2
(1 + δρ)
= −#
,
∂ρ
ρ + δρ2
∂λ1
= 1,
∂u
∂λ2
= 1.
∂u
Gli invarianti di Riemann devono quindi soddisfare le equazioni
4
4
"
"
1 + δρ
1 + δρ
λ+ = u + dρ
,
λ− = u − dρ
ρ
ρ
(A.27)
(A.28)
(A.29)
Procedendo con il calcolo degli integrali si distinguono i casi i cui α è maggiore
o minore si zero, ovvero η = ±1.
Se η = −1 gli invarianti di Riemann sono
.
,#
1
1
λ± = u ±
ρ − αρ2 − √ arctan
−1 ,
(A.30)
αρ
α
invece, se η = +1, si ha invece
λ± = u ±
,#
,
#
1
ρ − αρ2 + √ ln 1 + 2αρ + αρ + α2 ρ2 ,
2 α
(A.31)
Nel limite α → 0, l’Eq.(4.1) diventa la NLS defocusing e gli invarianti di
√
Riemann (A.30) e (A.31) si riducono a λ± = u ± 2 ρ, risultati già noti in
letteratura [10] .
Appendice B
Calcolo di onde cnoidali per
la NLS Quintica
Partendo dall’equazione (4.48), ponendoci nel caso η = −1, si vuole ricavare
la forma dei solitoni dark e antidark. Raccogliendo il termine 4/3α, si può
riscrivere la (4.48) come:
$ dρ %2
dθ
con Q(ρ)
Q(ρ) =
*
− ρ4 +
=
4
αQ(ρ)
3
+
3 3 6
6
3
ρ + k0 ρ2 + E0 ρ + C02 .
α
α
α
α
(B.1)
(B.2)
Supponendo quattro zeri distinti della Q(ρ) ordinati nel seguente modo e1 <
e2 < e3 < e4 , la (B.1) si può riscrivere
$ dρ %2
dθ
=
87
87
87
8
4α 7
e4 − ρ ρ − e3 ρ − e2 ρ − e1 ,
3
(B.3)
i quattro zeri sono legati alle costanti di integrazione tramite:
3 2
C ,
α 0
6
e1 e2 e3 + e1 e2 e4 + e1 e3 e4 + e2 e3 e4 = − E0 ,
α
6
e1 e2 + e1 e3 + e1 e4 + e3 e4 + e2 e4 = k0 ,
α
3
e1 + e2 + e3 + e4 = .
α
e1 e2 e3 e4 =
Procedendo con l’integrazione per quadratura si ricava
.
"
"
3
dρ
#
dθ =
.
4α
Q(ρ)
(B.4a)
(B.4b)
(B.4c)
(B.4d)
(B.5)
Si devono ora legare gli estremi di integrazione della (B.5) al potenziale Q(ρ).
Graficando il potenziale per quattro valori generici di e4 , e3 , e2 e e1 (B.1) si
nota che ci sono due zone in cui Q(ρ) > 0 e quindi due possibili soluzioni delle
67
68
0.2
0.5
1.0
1.5
2.0
!0.2
!0.4
!0.6
Figura B.1: Andamento del potenziale Q(ρ) per quattro zeri generici.
equazioni del moto. Considerando dapprima la zona in cui Q(ρ) > 0 per valori
di ρ compresi tra e1 < ρ < e2 , si ha
.
" θ
" ρ
3
dρ
#
dθ =
.
(B.6)
4α
Q(ρ)
θ0
e1
La risoluzione dell’integrale (B.6) si ottiene operando la sostituzione [9]
sin2 φ =
e4 − e2 ρ − e1
.
e2 − e1 e4 − ρ
(B.7)
Invertendo tale relazione si ottiene immediatamente un fattore del potenziale
Q(ρ)
(ρ − e1 ) = k̃(e4 − ρ).
(B.8)
dove per semplificare si definisce
K̃ ≡
e2 − e1
sin2 φ,
e4 − e2
(B.9)
Tramite alcune manipolazioni della (B.7), si possono ottenere gli altri fattori
del potenziale
7
87
e1 − e3 8
(ρ − e3 ) = e4 − ρ k̃ +
,
e4 − ρ
7
87
e1 − e2 8
(ρ − e2 ) = e4 − ρ k̃ +
,
e4 − ρ
e1 − e4
(e4 − ρ) =
.
1 + k̃
(B.10a)
(B.10b)
(B.10c)
Sostituendo le (B.10) nel potenziale G(ρ), si ottiene:
*
+2 4*
+*
+*
+
9
e4 − e1
e4 − e2
e2 − e1
e3 − e1
G(ρ) =
cos φ sin φ (1 − m1 sin2 φ)
e2 − e1
e1 − e4
e1 − e4
1 + k̃
(B.11)
69
dove si è definito il modulo m1 come
m1 =
(e4 − e3 )(e2 − e1 )
.
(e3 − e1 )(e4 − e2 )
(B.12)
Invertendo la relazione (B.7), si ricava l’espressione di ρ in funzione in φ
8 2
7
1
e1 + e4 ee24 −e
−e2 sin φ
8 2
7
,
(B.13)
ρ=
1
1 + ee24 −e
−e2 sin φ
differenziando si ricava
dρ =
*
e2 − e1
e4 − e2
+
*
(e4 − e1 ) sin 2φ
+2 .
7 e −e 8 2
1 + e24 −e12 sin x
Sostituendo (B.14), (B.11) nella (B.5) si ottiene
.
3
2F (φ, m1 )
#
θ − θ0 =
.
4α (e3 − e1 )(e4 − e2 )
Definendo
(B.14)
(B.15)
.
α(e3 − e1 )(e4 − e2 )
(B.16)
3
si ricava dalla (B.13) l’espressione dell’onda cnoidale relativa all’integrale tra gli
zeri e1 e e2
7
8 2
1
e1 + e4 ee24 −e
−e2 sn (A(θ − θ0 ), m1 )
8
ρ=
.
(B.17)
1
2
1 + ( ee24 −e
−e2 sn (A(θ − θ0 ), m1 )
A≡
L’integrale nella seconda area, compresa tra e3 e e4 , equivale alla risoluzione
dell’equazione
.
" θ
" e3
3
dρ
#
dθ =
.
(B.18)
4α
Q(ρ)
θ0
ρ
ottenuta, in modo analogo, tramite la sostituzione
sin2 φ =
(e3 − e1 )(e4 − ρ)
,
(e4 − e3 )(ρ − e1 )
(B.19)
ottenendo l’onda cnoidale
dove m2 è pari a
7
8 2
3
e4 + e1 ee43 −e
−e1 sn (A(θ − θ0 ), m2 )
7
8
ρ=
.
3
2
1 + ee43 −e
−e1 sn (A(θ − θ0 ), m2 )
(B.20)
(e4 − e3 )(e4 − e2 )
,
(e2 − e1 )(e3 − e1 )
(B.21)
e4 + re1 tanh2 A(θ − θ0 )
1 + r tanh2 A(θ − θ0 )
(B.22)
m2 =
Il limite di solitone è legato all’imposizione che m → 1, che comporta e3 = e2 =
ρ0 nelle onde cnoidali. Si ottiene le espressione del solitone dark
ρ1 =
70
e del solitone antidark
ρ2 =
e1 + r" e4 tanh2 A(θ − θ0 )
,
1 + r" tanh2 A(θ − θ0 )
(B.23)
dove r e r" sono rispettivamente
r=
e2 − e1
e4 − e2
r" =
e4 − e3
.
e3 − e1
(B.24)
Bibliografia
[1] G.B. Whitham, Proc . Roy. Soc. A 283, 238-261 (1965).
[2] G.B. Whitham, Linear and Nonlinear Waves, (Wiley-Interscience, 1974).
[3] A.V. Gurevich e L.P. Pitaevskii, Sov. Phys. JETP 38, 291 (1974).
[4] F. Calogero e A. De Gasperis, Spectral Transform and Solitons I (NorthHolland Publishing Company, 1982)
[5] H. Flaschka, M. G. Forest, and D. W. McLaughlin, Multiphase averaging
and the inverse spectral solution of the Korteweg-de Vries equation, Comm.
Pure Appl. Math., 33:739-784 (1980).
[6] G. A. El, V. V. Geogjaev, A. V. Gurevich, e A. L. Krylov, Physica D 87,
186 (1995)
[7] G.A. El, Chaos, 15, 037103 (2005).
[8] A.M. Kamchatnov, Nonlinear Periodic Waves and Their Modulations
(Singapore, World Scientific, 2000).
[9] P.F. Byrd e M.D. Friedman, Handbook of Elliptic Integrals for Engineers
and Physicists (Springer-Verlag, 1954).
[10] A. V. Gurevich and A. L. Krylov, JETP 92, 1684 (1987).
[11] D. J. Kaup, Prog. Theor. Phys. 54, 396?408 (1976)
[12] R. K. Dodd, J. C. Eilbeck, J. D. Gibbon, e H. C. Morris, Solitons and
Nonlinear Wave Equations (Academic Press, New York, 1982).
[13] C. De Angelis, IEEE J. Quantum Electron. 30 818 (1994)
[14] F. Kh Abdullaev, A. Gammal, L Tomio e T. Frederico, Phys. Rev. A 63
043604 (2001)
[15] Yu. S. Kivshar e G. P. Agrawal, Optical Solitons: From Fibers to Photonic
Crystals (San Diego, Academic, 2003)
[16] S. Gantz e J. Herrmann, J. Opt. Soc. Am. B, 8, 2296 (1991)
[17] I. V. Barashenkov, Phys. Rev. Lett. 77, 1193 (1996).
[18] W. Wan, S. Jia, e J. W. Fleischer, Nature Phys.03, 46 (2007).
71
Bistability and instability of dark-antidark solitons
in the cubic-quintic nonlinear Schrödinger equation
M. Crosta
Dept. of Physics, Sapienza University of Rome, I-00185, Rome, Italy
A. Fratalocchi
PRIMALIGHT, Faculty of Electrical Engineering; Applied Mathematics and Computational Science,
King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia∗
S. Trillo
Dipartimento di Ingegneria, Università di Ferrara, Via Saragat 1, 44122 Ferrara, Italy†
(Dated: August 11, 2011)
∗
†
O
TI
TR
IS
D
R
N
O
T
FO
RE
V
IE
W
Spatial solitons are important for their capability to
beat diffraction and their potential for engineering a variety of optical reconfigurable structures including (and not
limited to) couplers, deflectors and logic gates. In Kerr
media where the paraxial propagation is described by
the nonlinear Schrödinger (NLS) equation, only one soliton solution exists once the parameters (i.e., nonlinearity
and peak intensity or width) are fixed and that solution
is stable. However, more general nonlinear responses can
result into bistability of solitons (strictly speaking solitary waves) and/or their instability against the growth
of weak perturbations.
In this paper we are interested to investigate such features for solitons sitting on a finite background (i.e., darklike) in the context of the Cubic-Quintic NLS (CQNLS)
with a defocusing cubic and focusing quintic nonlinear
response. The importance of such model lies in the fact
that it constitutes the simplest model for a defocusing
saturable Kerr effect [1, 2], whose parameters can be effectively measured in a relatively simple way by two-wave
coupling or Z-scan [3, 4].
As far as bistability is concerned, the case of dark
solitons has been investigated with reference to various
model including the CQNLS [2, 5–9]. In particular Hermann has shown that dark solitons exhibit bistability of
the second kind, i.e. characterized by solutions possessing the same full-width-half-maximum (FWHM), albeit
possessing different amplitudes (and generally invariants
of motion). This type of bistability was introduced in
Refs. [10, 11] by Gatz and Hermann, to distinguish it
from the earlier definition [12–14] which implies the existence of different solutions possessing the same value
U
of one invariant of motion (e.g., the power) for different
values of the internal parameter, typically the nonlinear
propagation constant β. The analysis carried out by Herrmann, however, is limited to stationary solitons, while a
full family of moving dark solitons can exist. The analysis is further complicated by the fact that the CQNLS in
the regime considered here is known to possess coexisting
antidark (or bright on pedestal) solutions [9, 15]. The full
family of dark and antidark solitons, once parametrized
by the velocity, exhibits intriguing features which have
been overlooked, and which we discuss below. Furthermore, whether the full family of the dark-antidark moving pairs is stable or not, and which are the instability
mechanisms is still an open problem. Our systematic investigation of these problems provides two answers: (i) it
shows that the criterium demonstrated by Barashenkov
[16] for dark solitons provides the correct exhaustive answer to the stability problem also for antidark solitons;
(ii) it clarifies that the instabilities are not of the collapsing type as conjectured in the previous literature (collapse can in principle occur even in 1+1D with nonlinear
mechanism such as the cubic-quintic nonlinearity), rather
involving new scenarios in proper regions of the parameter space.
IB
INTRODUCTION
CO
I.
PY
N
We characterize the full family of soliton solutions sitting over a background plane wave and ruled
by the cubic-quintic nonlinear Schrödinger equation in the regime where a quintic focusing term
represents a saturation of the cubic defocusing nonlinearity. We discuss existence and properties
of solitons in terms of catastrophe theory and fully characterize bistability and instabilities of the
dark-antidark pairs, revealing new mechanisms of decay of antidark solitons.
http:www.primalight.org; [email protected]
[email protected]
Besides being important per se, the knowledge of the
dynamics of the whole soliton family of the CQNLS is
also important in view of recent studies which extend
the investigation of competing nonlinearities to the nonparaxial [17] and nonlocal [18, 19] regimes. Moreover, the
full characterization of the soliton solutions and their instabilities constitute the starting ground for describing
the feature of dispersive shock waves (DSW, involving
multiple solitons in the weakly dispersive regime) [20],
an active area of research where succesful experiments
have been recently performed in non-Kerr media under
different excitation conditions [21–23]. In this respect,
here we provide the first prediction of a dispersive shock
wave produced directly by the decay of a solitary wave
2
∂u 1 ∂ 2 u
α 4
2
+
− |u| u + |u| u = 0,
2
∂z
2 ∂x
2
(1)
which describes a saturable Kerr-like nonlinearity
through its truncated expansion at second-order in the
2
normalized intensity |u| , with α being an external free
parameter that weights the quintic nonlinear response
(the smaller α, the weaker the saturation effect). Solitons of such system have been reported before. However, we reformulate the full problem from the beginning, giving novel analytical formulas which prove convenient for the purpose of our analysis. Solitons correspond to translationally
invariant solutions of the the
!
form u(z, x) = ρ(θ) exp [iφ(θ) + iβz], where θ = x − vz
and β = g(ρ0 ) is the nonlinear phase shift experienced by
a plane-wave background with intensity ρ0 ≡ |u0 |2 , in the
medium where the nonlinear refractive index varies with
2
intensity ρ = |u| according to the law g(ρ) = −ρ+αρ2 /2.
These solutions depend on two internal parameters which
we choose, in analogy to general dark soliton solutions
of the defocusing NLS equation, as ρ0 (intensity background, which fixes also β) and v (soliton velocity, which
fixes also the darkness or brightness of the soliton). Note
that, here, the quintic term prevents the simple rescaling
to ρ0 = 1 without rescaling α, and the velocity complicates further the scenario, so we keep the three parameters free. The modulus ρ plays the role of equivalent
”position”, and obeys the standard Hamiltonian dynamics with ”momentum” p = ρ̇ ≡ dρ/dθ,
∂H
∂H
p2
, ρ̇ =
; H=
+ V (ρ),
∂ρ
∂p
2
(2)
#
"
2
α 3
c
V = 2ρ (ρ − ρ30 ) − (ρ2 − ρ20 ) + 2k0 (ρ − ρ0 ) − 0 .
3
ρ0
ṗ = −
2
Here c0 = vρ0 and k0 = v2 + ρ0 − α2 ρ20 . Once ρ(θ) is
obtained by solving Eqs. (2), the phase profile φ(θ) can
be found by integrating the following equation:
$
%
ρ0
φ̇ = v 1 −
.
(3)
ρ
A.
Soliton solutions
Soliton solutions sitting on the plane-wave background
ρ = ρ0 correspond to homoclinic separatrix trajectories
of Eqs. (2), characterized by the energy level H = E with
E = V (ρ0 ) = −2c20 . Such separatrices emanate from the
saddle point (ρ, p) = (ρ0 , 0) of the Hamiltonian H(ρ, p).
ρd (θ) =
ρm + rρa tanh2 [w(θ − θ0 )]
,
1 + r tanh2 [w(θ − θ0 )]
(4)
where r = (ρ0 − ρm )/(ρa − ρ0 ) and w =
!
α(ρa − ρ0 )(ρ0 − ρm )/3 is the inverse soliton width,
while the minimum intensity (dip) is given by the root
ρ = ρm . Similarly for anti-dark solitons, we obtain
ρa + 1r ρm tanh2 [w(θ − θ0 )]
ρa (θ) =
1+
1
r
tanh2 [w(θ − θ0 )]
,
(5)
where r and w are the same as for dark solitons.
2
1
ρ
1
ρm
0
ρa
2
i
|u| , φ (rad)/π
We start from the CQNLS equation:
0
0.5
(a)
−1
0
1
2
3
(c)
0
4
−5
ρ
5
5
0
x
5
(d)
0
−5
0
x
3
(b)
|u|2, φ (rad)
DARK-ANTIDARK SOLUTIONS
V(ρ)−E
II.
For α #= 0, the potential V (ρ) − E is a double well corresponding to a double-loop separatrix as shown in Fig.
1(a,b). Therefore one has, in general, a coexisting pair of
dark and antidark solitons that corresponds to the motion along the left well ρm ≤ ρ ≤ ρ0 (dark solitons), and
the right well ρ0 ≤ ρ ≤ ρa (antidark solitons), respectively. Here ρm ≤ ρ0 ≤ ρa are the roots of V (ρ) − E
(explicit expressions of ρm , ρa are reported in Appendix
A). In terms of such roots, we derive (see Appendix A)
the following explicit solutions for dark solitons:
p
of the CQNLS model.
0
1
2
ρ
3
4
2
1
0
−5
FIG. 1. (Color online) Dark-antidark pair sitting on the unit
background ρ0 = 1, for α = 0.6, v = 0.3: (a) potential V (ρ) −
E vs. ρ; (b) phase-space picture (contour lines of H); (c-d)
Relative intensity (thick black solid line) and phase (thin red
solid line) profiles of dark (c) and antidark (d) solitons.
The expressions in Eqs. (4)-(5) allow us to obtain,
by integrating Eq. (3), the nonlinear phase associated to
the two soliton families. We obtain for dark and antidark
solitons, respectively
$&
%
v
rρa
tanh[w(θ − θ0 )] + φ0 , (6)
φd = − s tan−1
w
ρm
$&
%
v
ρm
φa = s tan−1
tanh[w(θ − θ0 )] + φ0 , (7)
w
rρa
'
(ρa −ρ0 )(ρ0 −ρm )
. From Eqs. (6-7) one can
where s ≡
ρa ρm
easily calculate the phase jump ∆Φ = Φ(+∞) − Φ(−∞)
across the soliton.
3
bounded by the curve (so called bifurcation set [24])
! a "3 # b $2
+
= 0,
(8)
3
2
ρ =1, v=0.5
0
−0.5
ρ0=3, v=0.5
ρ0=1, v=0.7
a
−1
ρ =3, v=0
0
−1.5
−2
ρ =1, v=0
0
−2.5
−1
ρ0=4, v=0
−0.5
0
b
0.5
1
FIG. 2. (Color online) Cusp catastrophe picture for darkantidark soliton pairs. Each curve shows the evolution of
the parameter a and b of the normal form potential V (y) =
y 4 /4 + ay 2 /2 + by, calculated for a dark-antidark soliton pair
with fixed internal parameters v and ρ0 , and α changing from
zero up to its critical value αc , where all the curves arrive
tangentially on the cusp curve (Eq. (8), red solid line) that
bounds the soliton existence domain.
2
(a)
0
0
−2
−2
−4
−4
(e)
0
0
1
2
3
4
4
0
1
2
3
4
4
(b)
p
0.1
(c)
−0.1
5
2
0
0.5
1
(d)
Taking fixed internal parameters ρ0 and v, while changing α continuosly, makes the control parameters a and
b calculated for the soliton to span a smooth curve in
the control parameter plane (a, b), until at α = αc , they
hit (arriving tangentially) the boundary set by the cusp
curve [Eq. (8)]. Different values of ρ0 and v result into
different control curves, as displayed in Fig. 2. We point
out that a similar behavior occurs by varying ρ0 or v,
keeping the other two parameters fixed. In particular,
in the latter case, the existence domain turns
%out to be
−vc ≤ v ≤ vc , with the cut-off velocity vc = ρ0 − αρ20 ,
obtained by expressing Eq. (9) in terms of v = vc .
1.5
1
(f)
2
2
1
0
0
0
−2
−2
−1
dip intensity
V(ρ)−E
2
shown in Fig. 2. Such curve marks the values where the
critical points (∂y V = 0) of the potential become doubly
degenerate (∂y2 V = 0), and exhibits the characteristic
shape of a cusp in the origin (three-fold degenerate point,
∂y3 V = 0). In terms of the original parameter ρ0 , v, α, the
existence condition requires α ≤ αc , with the following
critical value of the quintic coefficient αc :
#
$
v2
1
1−
.
(9)
αc =
ρ0
ρ0
(a)
0.6
0.3
peak intensity
0
0.8
0.5
α=0
2
10
0.1
1
10
0
1
2
ρ
3
4
−4
0
1
2
ρ
3
4
5
−2
0
0.5
1
1.5
ρ
FIG. 3. (Color online) Potential and phase-space picture for
v = 0, ρ0 = 1, and different values of α: (a-b) α = 0 (ideal
Kerr case); (c-d) α = 0.5; (e-f) α = 0.9 (close to critical value
αc = 1, where the right branch of the separatrix becomes
vanishingly small).
Interestingly, the existence domain of the soliton pairs
can be described with the aid of catastrophe theory
[24], already applied to characterize dark-antidark soliton pairs in a different context (gap soliton theory [25]).
In fact, the quartic potential V in Eqs. (2) belongs
to the A+ family, and gives rise to the so called cusp
catastrophe. According to this picture the potential
V (ρ) in Eq. (2), which is of the general form V (ρ) =
c4 ρ4 + c3 ρ3 + c2 ρ2 + c1 ρ, can be cast into the canonical
form [24] V (y) = y 4 /4 + ay 2 /2 + by, by means of the
change of variable ρ = (4c4 )−1/4 y − c3 /(4c4 ). Then, in
the control parameter plane (a, b) (explicit expressions
of a and b as a function of α, ρ0 , v are cumbersome but
can be easily derived), solitons exist in the inner region
0
0
0
0.5
velocity v
1
0.3
0.6
0.9
−4
α=0.01
(b)
10
0
0.9
0.5
velocity v
1
FIG. 4. (Color online) (a) dip intensity (darkness) of dark
solitons and (b) maximum intensity of antidark solitons as a
function of velocity v for different values of α. Here the background is ρ0 = 1. The dashed curves represent the existence
threshold set by Eq. (9). Note the vertical log scale in (b).
As an example, we show in Fig. 3 how the typical
phase plane (potential) changes when α is varied between zero and the critical value αc . In this case we
choose ρ0 = 1 and still solitons, viz. v = 0, yielding
ρm = 0 which means that the dark soliton is black regardless of the value of α, while the antidark is characterized by a peak intensity ρa = 3/α − 2ρ0 . In the
limit α = 0 shown in Fig. 3(a,b), which represents the
ideal Kerr case, the potential is cubic, and the separatrix has only one branch corresponding to the well-known
black soliton solution (ρ = tanh2 (x)) of the NLS equation, whereas for ρ > 1 the motion is unbounded and
no coexisting antidark solutions do exist. In fact the unbounded motion for ρ > 1 can be thought of as the motion in the right well that, however, becomes infinitely
4
deep [V (θ = ∞) → −∞)] and wide (since ρa → ∞).
Viceversa, as α grows from zero, the behavior of the potential at θ = ∞ is inverted, and the right well becomes
finite, allowing for a eight-shaped separatrix corresponding to the dark-antidark pair [see Fig. 3(c-d)]. For small
values of α the antidark soliton has high peak intensity ρa
above the background ρ0 , which, however, decreases as
the saturation parameter α increases. For α approaching
its critical value αc antidark solitons become shallow [see
Fig. 3(e-f)], until they reduce to the plane wave exactly
at α = αc , where ρa → ρ0 . The behavior of dark-antidark
FIG. 5. (Color online) Color level plot of FWHM [θF W HM in
Eq. (10)] in the parameter plane (α, ρ0 ) for v = 0.1, dark (a)
and antidark (b) solitons. (c,d) dark solitons with v = 0.4 (c),
and v = 0.8 (d). The dashed lines give the critical condition
αc (ρ0 , v).
solitons with velocity exhibits intriguing features, which
can be gathered by plotting the minimum (dip) intensity
ρm of dark solitons (the larger ρm , the lower the darkness) and the peak intensity ρa of antidark solitons (the
higher ρa , the brighter the antidark) versus v at constant ρ0 , for different values of α, as displayed in Fig.
4. For non-zero but small velocities the picture remains
qualitatively unchanged with respect to the case v = 0,
in the sense that the peak intensity ρa of antidark solitons decreases continuously from infinity (at α = 0) to
ρa = ρ0 at the critical value αc , such that they become
infinitely shallow (i.e. they reduce to a pure plane wave).
In this case, however, the dip intensity of dark solitons
is no longer zero, i.e. they becomes gray solitons with
darkness ρ0 − ρm . The darkness decreases for growing
velocities !
v up to a minimum value at the bound velocity vc = ρ0 − αρ20 (obtained by solving Eq. (9) with
respect to v for fixed α). This is clearly shown in Fig.
4, where we summarize the result for a fixed background
ρ0 = 1. Note from Fig. 4(a) that Kerr (NLS) dark solitons (α = 0) are always darker than the corresponding
CQNLS solitons of the same velocity (the curve ρm (v)
for α = 0 is always below the other curves relative to
α $= 0), and have also a larger phase jump ∆Φ than their
CQNLS counterparts.
Interestingly, however, as the velocity grows large
enough (above v = 0.5 in Fig. 4), it turns out that
dark solitons can become infinitely shallow at the cut-off
condition for their existence (ρa → ρ0 and hence darkness tends to zero). Conversely, under the same conditions, antidark solitons cease to become infinitely shallow, rather reaching a finite minimum peak intensity at
the cut-off condition for their existence [see Fig. 4(b)].
From Fig. 4(b) it is also clear that, for small α the brightness of antidark solitons is nearly independent on the
velocity. This change of behavior at the cut-off condition for the existence (from infinitely shallow antidark to
infinitely shallow dark solitons) depends the background
ρ0 . It can be shown to occur at the value of α0 = 3/(4ρ0 ),
√
in correspondence of the velocity v0 = ρ0 /2 (in Fig. 4,
α0 = 0.75 and v0 = 0.5).
The coexistence of dark and antidark solitons constitute bistable mechanism such that two different solutions exist with same parameters (ρ0 and v) and different renormalized invariants M, H, P (see Appendix A
for their definition). However these soliton families can
be bistable also according to the definition by Gatz and
Herrmann, i.e. for fixed α different solutions of the same
width can exist although they sit on a different background ρ0 . This type of bistability was investigated for
still (v = 0) dark solitons [7]. In order to generalize this
result to the full family (any v, and antidark case), we
have calculated the FWHM as
θF W HM =
!
2
tanh−1 f (ρ0 , ρm , ρa ),
w
(10)
a −ρ0
for dark solitons
where f (ρ0 , ρm , ρa ) = 2ρaρ−ρ
0 −ρm
(FWHM taken at half the intensity between ρ0 and the
ρ0 −ρm
dip ρm ) and f (ρ0 , ρm , ρa ) = ρa +ρ
for antidark soli0 −2ρm
tons (FWHM at half the intensity between ρ0 and the
peak ρa ), respectively. The results obtained by mapping
θF W HM in the plane (α, ρ0 ) are summarized in Fig. 5.
As shown, dark solitons exhibit bistability for sufficiently
large α, regardless of their velocity. However the range
of values of α where bistability occurs is greatly reduced
for large velocities [see Fig. 5(c)]. Conversely, as shown
in Fig. 5(d), antidark solitons are never bistable in the
same sense (no folding of the level curves of θF W HM is
ever observed at any velocity).
III.
INSTABILITY SCENARIOS
Having characterized the features of dark-antidark soliton pairs, we discuss their stability. We proceed by applying a known stability criterium according to which the
stability is related to the derivative of the invariant momentum M (see Appendix A) of the soliton against its
5
2
1
(b)
Hamiltonian H
momentum M
(a)
1.5
dark
1
0.5
0
0
vc
antidark
0.3
velocity v
0.8
antidark
0.6
dark
0.4
0.2
0
0
0.6
0.5
1
1.5
momentum M
2
FIG. 6. (Color online) (a) Renormalized momentum M vs.
soliton velocity v, and (b) Hamiltonian H(v) vs. M (v), for
dark and antidark solitons with fixed ρ0 = 1 and α = 0.7.
The red bullets mark the marginal condition ∂v M = 0 for
antidark solitons.
0.6
1
α0
(a) ρ0=1
coefficient α
coefficient α
1
0.8
αm
0.4
0.2
0
0
0.8
(b) ρ0=2
0.6
0.4
α0
slope. Therefore the whole family of dark solitons is stable in its existence domain, a conclusion that is fully supported by our numerical simulation of the propagation.
Viceversa instabilities take place for antidark solitons.
In particular, for small α, it turns out that they are always unstable in the whole domain of existence since
M (v) exhibits always positive slope. Conversely at sufficiently large α, the momentum M (v) changes its slope
near the cut-off value for existence vc . This is shown in
Fig. 6, where we compare the momentum M (v) for dark
and antidark solitons at α = 0.7 and ρ0 = 1. Note that
the change of slope of the momentum means that two
antidark solutions with same momentum and different
velocity exist. These two solutions differ by their Hamiltonian H (see Appendix for its expression). Indeed, as
shown in Fig. 6(b) the Hamiltonian H as a function of
M is folded, showing an upper (unstable) branch and a
lower (stable) branch.
αm
0.2
0.5
1
velocity v
1.5
0
0
0.5
1
velocity v
1.5
FIG. 7. (Color online) Stability maps for antidark solitons
in the parameter plane (v, α) for (a) ρ0 = 1; (a) ρ0 = 2.
The light (cyan) and dark (red) shaded domains correspond
to unstable and stable solutions, respectively. The dot-dashed
line sets the value α0 below which antidark solitons have finite
brightness at cut-off. The dynamics illustrated in Fig. 8 and
Figs. 9-10 are relative to the sampled values marked by bullets
and stars in (a), respectively.
velocity v. The marginal condition
∂M
= 0,
∂v
(11)
separates stable solutions (∂v M < 0) from unstable ones
(∂v M > 0). Such criterium, proved by Barashenkov [16]
for dark solitons of the generalized NLS equation is found
to account also for the instabilities of antidark solitons.
We emphasize, however, that such criterium accounts for
real eigenvalues crossing into the right-half plane (usually
bifurcating from zero eigenvalues associated with neutral
modes, i.e. symmetries of the model). Other instability mechanisms such as oscillatory instabilities resulting
from edge bifurcations leading to pairs of complex conjugate eigenvalues entering the right-half plane, need an
independent characterization. However, we have found
(numerically) no evidence for instabilities of such kind.
Moreover, the condition for the existence of solitons [Eq.
(9)] turns out to coincide with the region where the background plane-wave is modulationally stable. Therefore
the condition (11) turns out to be an exhaustive criterium
to assess the linear stability of all the solutions presented
so far.
The calculation of the function M (v) for dark soliton
family shows that such function has always a negative
FIG. 8. (Color online) Dynamics of antidark solitons with
ρ0 = 1, α = 0.7, exhibiting different behavior depending on
the initial velocity v [see bullets in Fig. 7(a)]: (a) v = 0.5,
stable propagation; (b) v = 0.4, decay into a nearby stable
antidark soliton; (c) v = 0, decay into a pair of antidark
shallow solitons with opposite velocities.
The velocity which gives the marginal stability (∂v M =
0) clearly depends on ρ0 and α. The results of Fig. 6
could be repeated for different values of ρ0 and α, and
summarized by drawing stability maps in the plane v, α
at constant ρ0 . This is shown in Fig. 7 for two different
values of the background ρ0 . As shown a relatively small
island of stability is found in the vicinity of the boundary
for existence, where they have small brightness. In par-
6
FIG. 9. (Color online) Decay of an unstable antidark solitons
(ρ0 = 1, v = 0, α = 0.1) into two symmetric dispersive shock
fans (trains of dark solitons): (a) color level plot of the intensity; (b) snapshot of intensity (black solid line) and phase
(red solid line) at z = 5. For comparison the input intensity
(renormalized to the maximum of the plot) is shown (solid
blue line).
FIG. 10. (Color online) As in Fig. 7 for v = 0.8.
ticular stability requires α > αm , where αm corresponds
to the vertex of the stability island [see Fig. 7]. In order
to test the validity of the marginal stability condition
we have made extensive simulations. Such simulations
confirm indeed that ∂v M = 0 gives the threshold for stability also for antidark solitons, and reveal basically two
different mechanisms of instability.
The first scenario is valid for relatively large α, such
that a stable range of velocities exists. As an example for illustration purpose we have chosen α = 0.7,
ρ0 = 1, which yields a range of stable velocities v =
(0.4644, 0.5477) [see Fig. 6(a)]. The dynamics for three
different values of v, marked by bullets in the map of Fig.
7(a), is displayed in Fig. 8. As shown in Fig. 8(a), for
v = 0.5 the antidark soliton propagates in stable fashion.
Conversely, when v is decreased just below the threshold
for instability [v = 0.4 in Fig. 8(b)], the initial soliton
decays into a different antidark soliton of higher velocity,
which is stable. However, for lower velocities, such that
the soliton lies far from the threshold for stability, the
decay into a pair of stable antidark soliton is observed,
a process which is accompanied by radiation. Note that
the instabilities are extremely long-range in both cases,
i.e. they become manifest only after tens of diffraction
lengths. As shown in Fig. 8(c), the two antidark solitons are symmetric in the limit v = 0 in such a way
that the initial zero momentum is conserved, otherwise
they appear to be asymmetric (case not shown). The
scenario illustrated in Fig. 8 holds also for different values of α, providing α > αm , i.e. a stable range of velocities exists. It is interesting to note that the value
ρ0 = 3/(4ρ0 ), which discriminates the fact that antidark
solitons can be or not infinitely shallow at cut-off, lies
in the region of stability for any ρ0 (see dot-dashed lines
in Fig. 7). Therefore antidark solitons with arbitrarily
small brightness that exist for α > α0 around cut-off,
are always stable, whereas antidark solitons which have
finite brightness at cut-off can propagate stably only for
α > αm .
The scenario discussed above change qualitatively
when α < αm , where the island of stability shrinks to
zero. In this case the decay instability towards other antidark solitons is forbidden because no stable solutions
exist. In this case, the only stable solutions are dark solitons and therefore the decay instability occurs towards
these solutions, even though their shape differ dramatically (indeed being ”opposite”) from the input antidark
shape. Importantly when α is small (weak saturation)
the antidark solitons possess large amplitude and large
power Pa and hence number of particles (see Appendix
for its definition). Under these conditions we have found
that the decay instability of the antidark leads to a DSW
(see Refs. [20–23] and references therein), i.e. an expanding region filled with fast oscillations which behave
asympotically as solitons. In the present case, starting
with a zero-velocity antidark soliton, the decay instability leads to two symmetric DSW fans, as displayed in
Fig. 9 for α = 0.1. In each of the two fans, the inner
edge is set by the darkest and slowest soliton, whereas
on the outer edge the fan is linked to the plane wave
through a train of solitons with progressively decreasing
darkness and increasing velocity, which become denser as
the background is approached. Although this behavior is
reminiscent of that ruled by the integrable NLS equation
(limit α = 0) in the semiclassical regime (i.e., nonlinearity much stronger than diffraction/dispersion) under
excitation of e.g. a gaussian on pedestal [21], it must
be emphasized that this is, to the best of our knowl-
7
edge, the first example where the dispersive shock waves
occurs directly through the decay of an unstable solitary wave of the system. Further characterization of the
shock fan need to develop the Whitham modulational
theory for the system, which is beyond the scope of this
paper. However, simple arguments based on the features
of solitons as shown in Fig. 4, let us predict that the
fan ruled by Eq. (1) is narrower as compared to the one
ruled by the integrable NLS equation under the same excitation. This is due to the fact that the soliton with
vanishingly small darkness, constituting the outer edge
of the fan, correspond to progressively reduced velocity
as α increases, as clearly shown in Fig. 4(a).
The decay scenario shown in Fig. 9 does not change
when starting with an antidark soliton with non-zero velocity (v, M != 0), except that the two shock fans become
asymmetric (the higher the velocity, the higher the asymmetry), in both the number of solitons and the velocity
of the darkest soliton (inner edge of the fan), as shown
in Fig. 10. This asymmetry is largely expected on the
basis that a symmetric configuration would have M = 0
and hence would lead to violation of momentum conservation. Moreover the asymmetric development of the
shock is analogous to DSW generated in the integrable
NLS limit when starting from a gray beam with non-zero
velocity (momentum) [23].
Finally we point out that, for larger α, yet with α < αc ,
the scenario remains qualitatively unchanged though the
number of dark solitons generated in the decay of the
antidark soliton decreases.
Though the aim of this paper was the full characterization of solitary waves sitting on a finite background, before concluding, we point out that such solitons can also coexist with bright solitary waves with
zero pedestal.
! Seeking for bright solitons of the form
u(x, z) =
ρ(x) exp(iβz), a simple calculation shows
that the peak!
intensity of these solutions turns out to
be ρa = 3(1 + 1 + 8βα/3)/(2α). The existence domain
of such solutions include arbitrarily small α, and hence
they coexist with dark-antidark pairs. Given the defocusing nature of the leading order (cubic) nonlinearity
this might appear surprising. However, it is not difficult
to understand that such bright solitons are sustained entirely by the quintic focusing nonlinearity since ρa diverges in the limit α → 0, where the quintic nonlinearity
vanishes indeed. Without deepening the study of the
bright case, already investigated in the framework of the
CQNLS equation, we expect them to be unstable at least
in the limit of small α, where they do not affect the decay
dynamics of antidark solitons or the stable dynamics of
dark solitons discussed above.
IV.
ity and instabilities have been studied for the full family
parametrized by the intensity of the background wave
and the velocity, and a new mechanism of instability involving the decay of an antidark soliton into a dispersive
shock wave has been characterized, and is expected to be
observable in saturable defocusing media.
V.
Analytical solutions can be worked out from the equation which follows directly from Eq. (2)
ρ̇ =
"
4α !
Q(ρ),
3
(12)
3
[E − V (ρ)]. The polynomial Q(ρ) can
where Q(ρ) = 4α
be expressed in terms of its ordered roots ρm , ρ0 (double
root), ρa , with ordering ρm ≤ ρ0 ≤ ρa , as Q(ρ) = (ρ −
ρa )(ρ − ρm )(ρ − ρ0 )2 . From Eq. (12) one obtains the
following quadrature integral,
"
# ρ(θ)
#
4α θ
dρ
!
=
dθ,
(13)
3 θ0
Q(ρ)
ρ(θ0 )
where ρ(θ0 ) = ρm or ρ(θ0 ) = ρa for dark or antidark
solitons, respectively. The integral (13) gives, upon inversion, the solutions given in the text. Explicitly, the
extremal roots ρm , ρa are expressed by:
!
3 − 2αρ0 − (3 − 2αρ0 )2 − 12v 2 α
,
ρm =
! 2α
3 − 2αρ0 + (3 − 2αρ0 )2 − 12v 2 α
ρa =
.
(14)
2α
Once the soliton solutions are known one can easily calculate the renormalized invariants (momentum, Hamiltonian) which are defined as follows
i
2
#
+∞
+∞
&
|ux |2
+
2
M=
H=
#
−∞
−∞
$
%
ρ0
(u∗x u − ux u∗ ) 1 −
dx,
ρ
(15)
'
(16)
#
ρ(x)
ρ0
[g(ρ0 ) − g(ρ)]dρ dx,
and the power or number of particles, Pa (antidark) and
Pd (dark):
CONCLUSION
Pa =
In summary we have discussed the main properties of solitary solutions with finite background of the
CQNLS equation with focusing quintic term. Bistabil-
APPENDIX A
#
+∞
−∞
|u|2 − ρ0 dx; Pd =
#
+∞
−∞
ρ0 − |u|2 dx(17)
The quantities M and H are those employed in the text
to assess the stability of soliton solutions.
8
[1] W. Krolikowski and B. Luther-Davies, Opt. Lett., 18,
188 (1993).
[2] Y. S. Kivshar and B. Luther-Davies, Phys. Rep., 298, 81
(1998).
[3] Y.-F. Chen, K. Beckwitt, F. Wise, B. G. Aitken, J. S.
Sanghera, and I. D. Aggarwal, J. Opt. Soc. Am B, 23,
347 (2005).
[4] Z. Q. Zhang, W. Q. He, C. M. Gu, W. Z. Shen, O. H.,
and Q. X. Guo, Appl. Phys. Lett., 91, 221902 (2007).
[5] I. V. Barashenkov and V. G. Makhankov, Phys. Lett. A,
128, 52 (1988).
[6] L. J. Mulder and R. H. Enns, IEEE J. Quantum Electron., 25, 2205 (1989).
[7] J. Herrmann, Optics Communications, 91, 337 (1992).
[8] F. G. Bass, V. V. Konotop, and S. A. Puzenko, Phys.
Rev. A, 46, 4185 (1992).
[9] S. Tanev and D. I. Pushkarov, Opt. Commun., 141, 322
(1997).
[10] S. Gatz and J. Herrmann, J. Opt. Soc. Am. B, 8, 2296
(1991).
[11] S. L. Eix and R. Enns, Phys. Rev. A, 47, 5009 (1993).
[12] A. E. Kaplan, Phys. Rev. Lett., 55, 1291 (1985).
[13] R. H. Enns and S. S. Rangnekar, Phys. Rev. Lett., 57,
778 (1986).
[14] R. H. Enns, S. S. Rangnekar, and A. E. Kaplan, Phys.
Rev. A, 35, 466 (1987).
[15] Y. S. Kivshar, V. V. Afanasjev, and A. W. Sneider, Opt.
Commun., 126, 348 (1996).
[16] I. V. Barashenkov, Phys. Rev. Lett., 77, 1193 (1996).
[17] J. M. Christian, G. S. McDonald, and P. ChamorroPosada, Phys. Rev. A, 81, 053831 (2010).
[18] E. N. Tsoy, Phys. Rev. A, 82, 063829 (2010).
[19] Z. Zhou, Y. Du, C. Hou, H. Tian, and Y. Wang, J. Opt.
Soc. Am. B, 28, 1583 (2011).
[20] M. A. Hoefer, M. J. Ablowitz, I. Coddington, E. A. Cornell, P. Engels, and V. Schweikhard, Phys. Rev. A, 74,
023623 (2006).
[21] W. Wan, S. Jia, and J. W. Fleischer, Nature Phys., 3,
46 (2007).
[22] C. Conti, A. Fratalocchi, M. Peccianti, G. Ruocco, and
S. Trillo, Phys. Rev. Lett., 102, 083902 (2009).
[23] A. Fratalocchi, A. Armaroli, and S. Trillo, Phys. Rev.
A, 83, 053846 (2011).
[24] R. Gilmore, Catastrophe theory for scientists and engineers (Dover Publications, Inc., New York, 1993).
[25] C. Conti and S. Trillo, Phys. Rev. E, 64, 036617 (2001).
Crossover dynamics of dispersive shocks in Bose-Einstein condensates characterized
by two and three-body interactions
1
M. Crosta1,3 , S. Trillo2 and A. Fratalocchi1∗
PRIMALIGHT, Faculty of Electrical Engineering; Applied Mathematics and Computational Science,
King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia
2
Dipartimento di Ingegneria, Università di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
3
Dept. of Physics, Sapienza University of Rome, I-00185, Rome, Italy
(Dated: November 10, 2011)
We show that the perturbative nonlinearity associated with three-atom interactions, competing
with standard two-body repulsive interactions, can change dramatically the evolution of 1D dispersive shock waves in a Bose-Einstein condensate. In particular, we prove the existence of a rich
crossover dynamics, ranging from the formation of multiple shocks regularized by coexisting trains
of dark and antidark matter waves, to 1D soliton collapse. For a given scattering length, all these
different regimes can be accessed by varying the number of atoms in the condensate.
PACS numbers: 03.75.Lm,03.75.Kk,42.65.Sf
Introduction. — Bose-Einstein Condensation (BEC)
has been successfully described in the mean-field limit
by the Gross-Pitaevskii equation (GPE), which embodies
quantum two-body interactions (s-wave scattering) into
its cubic nonlinear term [1]. Though not compulsory for
(ideal, i.e. interaction-less) condensation, such term sustains several, remakable and otherwise unobservable, coherent phenomena such as soliton formation [2], quantum
shocks [3, 4], collapse [5], nonclassical states [6], time reversal, quantum chaos and Anderson localization [7]. The
exceptionally huge tunability (in magnitude and sign) of
the scattering length a through Feshbach resonances [8],
have made the area of ultracold atoms a prolific ground to
investigate new quantum and nonlinear physics. Particular interest, in this context, was stirred by higher-order
few-body interactions [9] boosted also by the observation
of resonant Efimov states [10]. Most efforts in this area
are presently aimed at assessing the impact of few-body
recombination loss coefficients which limits the lifetime
of condensates [11]. However, almost unexplored is also
the effect of three-body elastic collisions, which naturally
arises as a quintic conservative (Hamiltonian) term in the
higher-order expansion of the GPE [9, 12, 13] whose impact clearly grows as higher densities are being reached
in experiments [14]. It is therefore of paramount importance to assess whether the three-body interactions can
lead to a clear signature in terms of qualitative new scenarios, especially if the latter becomes accessible in the
perturbative limit.
In this letter we address this question with reference to
the dynamics of dispersive shock waves (DSWs), which
form in a repulsive BEC when the kinetic spreading
regularizes the tendency driven by the nonlinearity to
form multivalued wave fronts. At variance with earlier
[3, 4, 15] and more recent studies on DSWs [16–19], which
are all based on the standard GPE, we investigate the
case characterized by an attractive three-body nonlinear correction to the repulsive s-wave scattering. This
regime is expected to be achievable by exploiting resonance tuning for bosons and turns out to be relevant
also for superfluid fermions [20]. Unlike the case of competing nonlinearities of the same sign (where three-body
contributions merely strengthen the nonlinearity at high
densities), we find that the presence of the quintic term of
different sign leads to a rich phase-diagram of the wavebreaking phenomenon, with novel regimes that uniquely
identify the contribution of the three-body interaction.
Remarkably, for a given scattering length, the crossover
between different dynamics can be totally controlled by
changing the number of atoms in the condensate. Due to
the ubiquitous nature of DSW, these findings are also of
fundamental interest in other areas of classical physics,
including nonlinear optics [21], electronic systems [22],
and granular chains [23].
Model and general shock scenario. – We start from
the following dimensionless cubic-quintic GPE, describing the free evolution in 1D (after trap is released along
x) of a BEC ruled by two and three-body interactions:
∂ψ !2 ∂ 2 ψ
+
− |ψ|2 ψ + α|ψ|4 ψ = 0,
(1)
∂t
2 ∂x2
!
having introduced dimensionless spatial x = X! mN ω/!
and temporal t = T N ω/! coordinates, with ! =
!/mωRT F " 1 (RT F is the Thomas-Fermi
radius), and a
!
rescaled wave function ψ" = Ψ 4π!a/mωN , being m the
single boson mass, N = |Ψ|2 dX the number of particles
and ω the transverse BEC trap frequency. The parameter
α/!ω = N g3 /|g2 |2 (with α ≥ 0) uniquely determines the
properties of the system in terms of the number of particles and the relative strength of the three-body term g3
[9] over the two-body term g2 = 4π!2 a/m. At variance
with the cubic GPE system (α = 0), which is completely
characterized by the scattering length a, this normalization clearly shows that the high-order BEC dynamics
lives in a two dimensional phase space encompassing both
a and N . This will permit to explore different scenarios
i!
2
FIG. 1. Color Online. Dispersionless evolution for ρ+
0 = 1.5:
(a-b) as obtained from Eqs. (5-7) with α = 0.15 (a), α = 0.3
(b); (c-d) numerical results from Eqs. (2) comparing the case
α = 0.4 (c) with the case of a standard GPE (α = 0) with
positive cubic nonlinearity (d).
even for a constant scattering length. We study the general shock dynamics in the dispersionless (or hydrodynamic) limit of Eq. (1), which can be obtained by !apply√
ing the Madelung transformation ψ = ρ exp[ !i dxu]
for # → 0:
∂ρ ∂(uρ)
+
= 0,
∂t
∂x
∂u
∂u ∂f
+u
+
= 0,
∂t
∂x ∂x
(2)
where f (ρ) = ρ − αρ2 , and ρ, u act as density and velocity of the fluid, respectively. The impact of the quintic term can be understood by studying the decay of
an initial jump in the density, which can be easily realized experimentally. This amounts to solve the Riemann
problem associated with the initial value u(x, 0) = 0,
+
ρ(x, 0) = ρ+
0 + (ρ− − ρ0 )Θ(x) [Θ(x) is the Heaviside
±
−
+
function], ρ0 (ρ0 > ρ0 ) being the initial constant densities across x = 0. We solve Eqs. (2) by deriving the
following Riemann invariants λ± = u ± R(ρ):
R(ρ) =
"
ρ(1 − 2αρ) +
cos−1 (1 − 4αρ)
√
,
2 2α
(3)
±
which transform Eqs. (2) into the diagonal form ∂λ
∂t +
"
∂λ±
V± ∂x = 0, where V± = u ± ρ(1 − 2αρ) are the Riemann characteristic eigenvelocities. As long as V± are
strictly real, i.e.
α < 1/(2ρ−
0 ),
(4)
the evolution of λ± can be given in terms of the selfsimilar variable ζ = x/t through the following simple
wave solutions of the system (characterized by one of the
two Riemann variables remaining constant)
"
R(ρ+
ρ(1 − 2αρ) = ζ,
0 ) − R(ρ) −
"
−
R(ρ) − R(ρ0 ) + ρ(1 − 2αρ) = ζ.
(5)
FIG. 2. Color Online. Global phase diagram of the shock
+
dynamics in the parameter plane (ρ0 ≡ ρ−
0 , α) [ρ0 = 1].
The two simple waves (5) are connected by an intermediate region with constant values ρ = ρi and u = ui .
The latter can be found by the matching equations
−
+
−
ui = [R(ρ+
0 ) − R(ρ0 )]/2 and 2R(ρi ) = R(ρ0 ) + R(ρ0 ).
The borders xi± of the intermediate region are then calculated by substituting the density ρi into Eqs. (5):
"
xi±
≡ ζi± = R(ρi ) − R(ρ−
ρi (1 − 2αρi ).
(6)
0)±
t
Conversely, the simple wave external edges xe± ≡ ζe± t
are calculated by matching Eqs. (5) with the asymptotic
(input) values ρ±
0 , thus obtaining
#
√
xe−
xe+
+
= − ρ+
= 1 − 2α.
(7)
0 (1 − 2αρ0 ),
t
t
Equations (5)-(7) predict a rich scenario for the shock
dynamics, as illustrated in Fig. 1a-b. Henceforth, with−
out loss of generality, we take ρ+
0 = 1 and ρ0 ≡ ρ0 . For
small values of α [see Fig. 1a], wave-breaking occurs only
for x > 0, where the velocity ζi+ of the point xi+ is larger
than ζe+ , leading to a multivalued region, whose dispersive regularization leads to a DSW. However, a more
complex dynamics is observed when α is increased [see
Fig. 1b]: in this case not only xi+ but also xe− moves
sufficiently fast and induces wave breaking for x < 0,
thus producing a double shock in the dynamics. The
crossover between these two regimes is given by the condition ζe− (α∗ ) = ζi− (α∗ ). It is worth noting that a
double shock have also been predicted in a BEC flowing through a penetrable barrier [17], though both the
underlying mechanism and the regularization are completely different in the latter case. A third regime settles
in when Eq. (4) is violated, resulting in imaginary eigenvelocities V± [Eqs. (2) are no longer hyperbolic]. In this
case, solutions (5)-(7) are meaningless and we resort to
numerical integration of Eqs. (2). Our analysis show
(Fig. 1c) that the system evolution is characterized by
the generation of two opposite velocities (Fig. 1c, u in
the inset), which compress the wavefront and generate
a cusp-like singularity at x ≈ 0− . This behavior originates from a three-body dominance over the two-body
3
FIG. 3. Color Online. (a) Velocity of DSW edges ζt , ζl , and
rarefaction edges ζe− , ζi− vs. α; (b-c) Color level plot of
density evolution ρ(x, t) from Eq. (1) for $ = 0.01 and (b)
α = 0, (c) α = 0.1; For comparison the edges of the DSW
xl , xt from Eq. (9), xe− , xi− of the rarefaction wave from
Eqs. (6)-(7) are reported. (d) Snapshots ρ(x, t = 8) from
(b-c) showing the DSW dark soliton train for α = 0, 0.1.
term, as demonstrated by comparing Fig. 1c to Fig. 1d,
which shows the evolution of the same input launched in
a standard GPE with positive nonlinearity [Eq. (1) with
α = 0 and x → ix]. The two wave-breaking are qualitatively identical, which can be understood by considering
that the GPE with attractive two-body nonlinearity exhibits complex eigenvelocities, in the same way as Eq.
(1) for α > 1/2ρ0 . Despite such similarity, the long-term
dynamics (regularization) in the presence of three-body
nonlinearities will be totally different, as discussed in the
following paragraphs. Figure 2 summarizes the results
of the dispersionless analysis in a wave-breaking phase
diagram. In particular, the two curves α∗ = α∗ (ρ0 ) and
α = 1/(2ρ0 ), divide the phase space into three distinct regions, each characterized by different wave breaking scenarios deepened below.
Real eigenvelocities: tuning the shock dynamics. –
We study the shock regularization for α ≤ α∗ by
exploiting Whitham theory of modulation, in the form
suitable for nonintegrable systems [24]. We summarize
only the outcome of this approach, while deferring the
(involved) mathematical details to a successive paper.
When α ≤ α∗ , the multivalued region for x > 0 (Fig.
1) is regularized by the formation of a single DSW such
that a modulated cnoidal wave (or dark soliton train)
appears within a shock fan limited by a leading xl and
a trailing xt edge. By matching the dispersionless limit
with the Whitham equations [24], we are able to find the
velocities of the edges ζl = xl /t and ζt = xt /t (replacing
the hydrodynamic estimates ζe+ , ζi+ ), in the form:
!
"
1
,
ζl = (1 − 2α) 2γl −
γl
#
ζt = R(ρi ) − R(1) + γt ρi (1 − 2αρi ),
(8)
FIG. 4. Color Online. (a) Snapshots of density ρ(x) obtained
from Eq. (1) with α = 0.3, $ = 0.01, and a input jump ρ0 =
1.5; (b)-(c) zoom of the boxed evolution in (a), showing the
opposite type of regularization taking place in the dynamics.
with γl = γ(1) and γt = γ(ρi ) arising from the solution
of the differential equation:
∂γ
(1 + γ)(1 + 2γ − 8γα)
=
,
∂ρ
2ρ(2γ + 1)(2αρ − 1)
(9)
integrated in ρ ∈ [1, ρi ], with the initial condition γ(1) =
1 for γl , and γ(ρi ) = 1 for γt . Figure 3a displays the
DSW edge velocities ζl , ζt [calculated by integrating Eq.
(9)] and those of the rarefaction wave ζe− , ζi− [from Eqs.
(6)-(7)] as a function of α. Increasing the three-body
contributions results into a nearly linear decrease of all
relevant velocities ζl , ζt , ζe− , ζi− , with a consequent reduction of the extension of both simple waves and the
shock fan. It is worth emphasizing that the two edges of
the DSW behave in a markedly different way: the leading
edge ζl is strongly influenced by α, while the variation of
ζt is much smaller. As a consequence, even when α is
increased by a small factor, the shock dynamics is significantly affected in terms of both shock fan and shock
overall angular direction. To verify these predictions, we
resort to numerical integration of Eq. (1). As shown in
Fig. 3b-c, a remarkable agreement is obtained between
theory and numerical simulations. The snapshot in Fig.
3d shows that, despite a small variation of α, a substantial reduction of the shock fan extension occurs (nearly
a factor of two for ∆α = 0.1) accompanied by a reduced
average velocity (center of mass closer to x = 0). In
summary, in the regime where a single shock is formed
for x > 0, all the features of the shock dynamics are significantly varied by slightly perturbing the strength of
the three-body interaction terms α.
Real eigenvelocities: multi-shock generation through
antidark solitons. – When α∗ < α < 1/2ρ0 , simulations
show the occurrence of a second DSW for x < 0 (Fig. 4),
in complete agreement with the hydrodynamic analysis.
In this case, however, the DSW for x < 0 differs substantially from the dynamics discussed above for x > 0.
The latter, in fact, is characterized by a cnoidal wave
composed by a train of dark-like oscillations (see Fig. 3),
4
FIG. 5. Color Online. Snapshots of density ρ vs. x from Eq.
(1) for α = 0.4 and # = 0.01: (a) shock regularization; (b)
antidark collapse instability.
owing to the condition ρ(xt ) > ρ(xl ) (see also Fig. 1a).
The DSW, more specifically, originates from a series of
oscillations that start from the leading edge xl (where the
shock matches the linear wave background at ρ = 1) and
culminate into a dark soliton at the trailing edge [Fig.
4c]. On the contrary, the DSW for x < 0 moves in the
opposite (backward) direction and has a larger density
on the leading edge, which yields ρ(xl1 ) > ρ(xt1 ), being
xl1 and xt1 the leading and trailing edge of the DSW at
x < 0, respectively. As a consequence, the modulated
wave train that regularizes the shock needs to be composed by bright entities, as shown in Fig. 4b. In this
case, the trailing edge is an antidark (bright on pedestal)
soliton solution of Eq. (1), which has been thoroughly
investigated recently [25]. Quite remarkably, although
these solitons are mostly unstable [25], the average velocity of the DSW for x < 0 is such that the soliton (and
hence the DSW) is totally stable on propagation.
Complex eigenvelocities: 1D BEC collapse. — When
α > 1/2ρ0 , the wave-breaking observed in the simulations of Eq. (1) reflects indeed the dominant character of three-body interactions (Fig. 5a-b). In the early
stage, a cusp-like singularity in x ∼ 0− is generated (Fig.
5a), as predicted by Eqs. (2). Contrary to the the case
of the standard attractive GPE, which shows a similar
early stage behavior [see Fig. 1d], the three-body dominance in Eq. (1) leads to a completely different long
term dynamics. In the quintic case, in fact, the singularity tends to evolve into antidark solitons, due to the
condition ρ(xl ) > ρ(xt ). Contrary to the previous case,
however, here the solitons velocity is nearly vanishing.
Zero-velocity antidark solitons of Eqs. (1) are always unstable [25], and they lead the BEC towards irreversible
collapse near x = 0 (Fig. 5b). We emphasize that such
a collapsing dynamics, which we have verified to occur
also for different inputs (e.g. Gaussian on pedestal), is
a unique features of Eq. (1). Collapse, in fact, cannot
be observed in the standard GPE (α = 0) that, owing to
its integrable nature, does not possess unstable solitons.
It is also worthwhile remarking the perturbative value of
the quintic term (α > 0.34) which induces such a dramatic and abrupt change in the dynamics.
In conclusion, we have shown that the perturbative nonlinearity arising from three-body elastic collisions can
dramatically alter the breaking scenario in repulsive
BEC. The generation of multiple DSWs involving antidark soliton trains, as well as 1D collapsing dynamics are
new peculiar behaviors that can be controlled by means
of the BEC number of atoms. These results are expected
not only to foster new perspectives in BEC physics, but
also to stimulate novel experiments in nonlinear optics,
where the quintic terms account for saturation of classical Kerr nonlinearities. We acknowledge funding from
PRIN 2009 project (No. 2009P3K72Z).
∗
[email protected]; www.primalight.org
[1] F. Dalfovo, S. Giorgini, L. V. Pitaevskii, and S. Stringari,
Rev. Mod. Phys. 71, 463 (1999).
[2] C. Becker et al., Nat. Phys. 4, 496 (2008).
[3] Z. Dutton, M. Budde, C. Slowe, L. V. Hau, Science 293,
663 (2001).
[4] M. A. Hoefer et al., Phys. Rev. A74, 023623 (2006).
[5] J. M. Gerton, D. Strekalov, I. Prodan, and R. G. Hulet,
Nature 408, 692 (2000).
[6] N. Bar-Gill, D. D. Bhaktavatsala Rao, and G. Kurizki,
Phys. Rev. Lett. 107, 010404 (2011).
[7] J. Martin, B. Georgeot, D. L. Shepelyansky, Phys. Rev.
Lett. 101, 074102 (2008).
[8] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev.
Mod. Phys. 82, 1225–1286 (2010).
[9] T. Köhler, Phys. Rev. Lett. 89, 210404 (2002). E.
Braaten and H.-W. Hammer, Phys. Rep. 428, 259 (2006).
[10] T. Kraemer et al. Nature 440, 315 (2006).
[11] P. A. Altin et al., Phys. Rev. A 84, 033632 (2011).
[12] F. Kh. Abdullaev et al., Phys. Rev. A 63, 043604 (2001);
F. Kh. Abdullaev and M. Salerno, Phys. Rev. A 72,
033617 (2005).
[13] W. B. Cardoso, A. T. Avelar, and D. Bazeia, Phys. Rev.
E 83, 036604 (2011).
[14] T. Shobu et al., Phys. Rev. A 84 033626 (2011).
[15] A. M. Kamchatnov et al., Phys. Rev. E 66, 036609
(2002); B. Damski, Phys. Rev. A 69, 043610 (2004);
V. M. Perez-Garcia et al., Phys. Rev. Lett. 92, 220403
(2004).
[16] J. J. Chang, P. Engels, and M. A. Hoefer, Phys. Rev.
Lett. 101, 170404 (2008).
[17] A. M. Leszczyszyn et al., Phys. Rev. A 79, 063608 (2009).
[18] R. Meppelink et al., Phys. Rev. A 80, 043606 (2009).
[19] G. Dekel, V. Farberovich, V. Fleurov, and A. Soffer,
Phys. Rev. A, 81, 063638 (2010).
[20] P. Pieri and G. C. Strinati, Phys. Rev. Lett. 91, 030401
(2003).
[21] J. E. Rothenberg and D. Grischkowsky, Phys. Rev. Lett.
62, 531 (1989); S. Jia, W. Wan, and J. W. Fleischer,
Phys. Rev. Lett. 99, 223901 (2007); C. Conti et al., Phys.
Rev. Lett. 102, 083902 (2009).
[22] E. Bettelheim, A. G. Abanov, and P. Weigmann, Phys.
Rev. Lett. 97, 246401 (2006).
[23] A. Molinari and C. Daraio, Phys. Rev. E 80, 056602
(2009).
[24] G. A. El, Chaos 15, 037103 (2005).
[25] M. Crosta, A. Fratalocchi, and S. Trillo, Phys. Rev. A
(to appear) http://arxiv.org/abs/1111.1872v1.