i UNIVERSITY OF BELGRADE FACULTY OF PHYSICS Nikola Z. Petrović EXACT SPATIOTEMPORAL TRAVELING AND SOLITARY WAVE SOLUTIONS FOR THE GENERALIZED NONLINEAR SCHRÖDINGER EQUATION Doctoral Dissertation Belgrade, 2013 ii УНИВЕРЗИТЕТ У БЕОГРАДУ ФИЗИЧКИ ФАКУЛТЕТ Никола З. Петровић ТАЧНА ТАЛАСНА И СОЛИТОНСКА РЕШЕЊА ГЕНЕРАЛИСАНЕ НЕЛИНЕАРНЕ ШРЕДИНГЕРОВЕ ЈЕДНАЧИНЕ докторска дисертација Београд, 2013 iii Информације о ментору и члановима комисије Ментор: др Миливој Белић, научни саветник Института за Физику, Универзитета у Београду, и редовни професор Тексас А&М универзитета у Дохи, Катар Чланови комисије: др Миливој Белић др Маја Бурић, редовни професор Физичког факултета Универзитета у Београду др Ђорђе Спасојевић, ванредни професор Физичког факултета Универзитета у Београду iv Захвалница Желео бих пре свега да се захвалим свом ментору проф. Миливоју Белићу за сву помоћ коју ми је пружио током свих ових година и без кога ова дисертација не би била написана. Такође бих се захвалио Министарству Просвете и Науке Србије под чијим пројектом ОИ 171-006 ја тренутно радим своја истраживања, као и Катарској Националној Фондацији за Истраживања која је асистирала моје истраживање под пројектима NPRP 25-6-7-2 и NPRP 09-462-1-074. Такође бих се захвалио групи проф. Белића: др Милану Петровићу, др Драгани Јовић, др Александри Стринић и другим члановима групе за сву помоћ. Посебна захвалност иде др Најдану Алексићу за изузетну помоћ у анализи стабилности решења добијених у дисертацији. Желео бих да се захвалим и студентима који су ме асистирали у истраживачком раду: Анасу Ал Бастамију, Хусеину Захредину и Моизу Бохри. Додатно бих се захвалио и члановима комисије за сво залагање у вези моје дисертације. Коначно, захвалио бих се својој супрузи Ташани за сву њену љубав и подршку, мојој сестри Катарини, мојим родитељима Зорану и Невенки, мојој баки Олги и деди Милану, мојој ташти Верици, мом зету Александру и свим мојим пријатељима и колегама. Ову тезу посвећујем мом сину Борису и свој мојој будућој деци. Они су моја мотивација и највећа инспирација. v Acknowledgement I'd like to first and foremost thank my mentor Prof. Milivoj Belić for all his assistance and help he has provided throughout all these years and without whom this Thesis would not have been made. I'd also like to thank the Ministry of Education and Science of Serbia under whose project OI 171-006 I am currently doing my research and the Qatar National Research Foundation which has assisted my research under the project NPRP 25-6-7-2 and NPRP 09-462-1-074. I'd also like to thank the group of Prof. Milivoj Belić: Dr. Milan Petrović, Dr. Dragana Jović, Dr. Aleksandra Strinić and other members of his team for all their assistance. Special thanks go to Dr. Najdan Aleksić for invaluable assistance regarding the stability analysis of the solutions obtained in this Thesis. I'd also like to thank the students I worked with and who have aided me in my research: Anas Al Bastami, Hussein Zahreddine and Moiz Bohra. In addition, I'd also like to thank the members of my committee for all their hard work involving my dissertation. Finally, I'd like to thank my wife Tašana for her unlimited love and support, my sister Katarina, my parents, Nevenka and Zoran, my grandparents, Olga and Milan, my mother-in-law Verica, my brother-in-law Aleksandar and all my friends and colleagues. I dedicate this Thesis to my son Boris and all of my future children. They are my motivation and my biggest inspiration. vi Апстракт Напредак у нелинеарној оптици умногоме зависи од наше способности да нађемо нова решења разних диференцијалних једначина које се природно јављају у системима где светлост интерагује са нелинеарном средином. Иако су рекреирање ових система кроз експеримент и компјутерска симулација система два најчешћа и плодотворна приступа, крајњи циљ остаје да се нађу егзактна решења ових система. Циљ ове тезе је да комбинује раније технике налажења егзактних решења диференцијалних једначина и примени их на нелинеарну Шредингерову диференцијалну једначину (НШДЈ). Конкретно, настао је недавно пробој у применама одређених техника експанзије у налажењу одређених егзактних решења НШДЈ. Упркос ограничењу у комбиновању решења због нелинеарности система и чињенице да не могу општа решења да се нађу, сама чињеница да можемо идентификовати нека егзактна решења је од великог значаја за област, посебно код евалуирања какве су појаве могуће у таквим системима. Ова теза ће се фокусирати на примену технике Ф-експанзије користећи се Јакобијевим елиптичним функцијама (ЈЕФ) да би се решиле разне форме НШДЈ са нелинеарношћу трећег степена. НШДЈ са нелинеарношћу трећег степена је од фундаменталне важности за област нелинеарне оптике јер описује путовање светлости кроз материјал са Керовом нелинеарношћу. Одређеним модификацијама технике Ф-експанзије можемо наћи егзактна решења за широку класу система. Системи које ја презентујем у тези имају одређен скуп заједничких особина. Све једначине имају једну лонгитудиналну промењиву, или просторну или временску, због параксијалне апроксимације, и до три трансферзалне димензије, такође или просротне или временске понаособ. Ако су све трансферзалне вариабле просторне vii онда суму њихових других извода множим са коефицијентом дифракције β, а ако је нека од варијабли темпорална, онда говорим о коефицијенту дифракције/дисперзије. Та два коефицијента (дифракција и дисперзија) могу да се нормализују у један до на знак. У случају аномалне дисперзије коефицијенти имају исти знак, а у случају нормалне дисперзије супротан знак. Осим ова два коефицијента редукована у један, имамо такође и коефицијент χ који одређује јачину нелинеарности трећег степена, и коефицијент γ који одређује добитак (за позитивно γ) или губитак сигнала у нашем систему. Први систем који ћу проучити у својој тези је стандардна НШДЈ са кубичном нелинеарношћу и у аномалној и у нормалној дисперзији. У случају нормалне дисперзије, циметрија између просторних и временске варијабле је разбијена. Добијам да се решења два система слично понашају. Ако се коефигијенти дифракције/дисперзије и нелинеарности модификују да буду синусоидални онда се може добити стабилни солитонски и путујући таласи. Први систем који ћу проучити у својој тези је стандардна НШДЈ са кубичном нелинеарношћу и у аномалној и у нормалној дисперзији. У случају нормалне дисперзије, циметрија између просторних и временске варијабле је разбијена. Добијам да се решења два система слично понашају. Ако се коефигијенти дифракције/дисперзије и нелинеарности модификују да буду синусоидални онда се може добити стабилни солитонски и путујући таласи. Затим ћу применити Ф-експанзију на системе са нелинеарношћу вишег степена, специјално на системе са кубично-квинтичном нелинеарношћу и са септичном нелинеарношћу. Презентујем, колико ми је познато, прва егзактна и стабилна решења за ове системе. Након тога ћу применити методу на једначину Грос-Питаевског (ГП). Примељујући методу на једначину ГП сусрећем се са Рикатијевом viii дигеренцијалном једначином која нема опште решење за многе системе. Ипак, презентујем широку класу решења за оне системе код којих се Рикатијева једначина може решити. За константну јачину екстерног поља и коефицијента дифракције решења се распадају, те је потребан екстерни добитак да би решења задржала интензитет, али за синусоидан облике две функције могу се добити стабилна решења. Додатна решења се могу добити за комплексније облике јачине екстерног поља и коефицијента дифракције. Такође у тези ћу анализирати ефекат линеарног потенцијала на НШДЈ са нелинеарношћу трећег степена. Додатак линеарног потенцијала је детаљно проучен и нова решења су нађена за овај систем. На крају ћу анализирати стабилност добијених решења. Добија се да су у већини случаја решења модулационо стабилна кад се користи менажирање дисперзије. Кључне речи Јакобијеве елиптичне функције нелинеарна Шредингерова једначина Грос- Питаевски Област Нелинеарна оптика Ужа област Примена математичке физике у нелинеарној оптици УДК број: 535.58 (043.3) ix Abstract The progress of the field of non-linear optics greatly depends on our ability to find solutions of various differential equations that naturally occur in the systems where light interacts with nonlinear media. Though re-creating the systems through experiment and performing computer simulations are the two most common and fruitful approaches, the ultimate goal remains to find exact solutions of these systems. The goal of this Thesis is to combine the work done in the field of finding exact solutions to certain classes of non-linear differential Schrödinger equations (NLSE). Most notably, there has been a breakthrough as of late in applying various expansion techniques in finding certain exact solutions to various NLSE. Despite the limitations of combining said solutions due to the non-linear nature of the solutions and the fact that not all solutions can be found using these techniques, the very fact that we can identify certain exact solutions is of tremendous importance to the field, especially when it comes to evaluating the kinds of functions and behavior that are possible within such systems. This Thesis will focus primarily on applying the F-expansion technique using the Jacobi elliptic functions (JEFs) to solve various forms of the NLSE with the cubic nonlinearity. The NLSE with a cubic nonlinearity is one of fundamental importance in the field of nonlinear optics because it describes the travelling of a light wave through a medium with a Kerr-like nonlinearity. Through certain modification of the technique we can find exact solutions in a very large class of systems. The systems I present in this Thesis will share a certain set of common properties. All of the equations I will tackle have a single longitudinal variable, either temporal or spatial, due to the application of the paraxial wave approximation, and up to three transverse dimensions, again both temporal and spatial. If all the transverse variables are spatial I x assign to the sum of their second derivatives a diffraction coefficient β whereas if one of them is temporal, I speak of the diffraction/dispersion coefficient. The two coefficients can be normalized into one, up to their sign. In the case of anomalous dispersion, the two coefficients have the same sign. In the case of normal dispersion, the two coefficients have the opposite signs. Apart from these terms which are present in the ordinary wave equation of linear optics, we also have the third order nonlinearity whose strength is determined by a parameter χ and we also have the term γ which describes the gain of loss of the signal inside our system. The first system I will tackle in the Thesis is the NLSE with a cubic nonlinearity, for both the anomalous and normal dispersion. In the case of the normal dispersion, the symmetry between the temporal and other transverse variables is broken, and a previously unknown ansatz had to be used. I obtained that the solutions for the two systems behave similarly. If one modifies the diffraction and nonlinearity sinusoidally, one can obtain stable solitary and traveling wave solutions. Then I will apply the F-expansion technique to obtain solutions for systems with higher- order nonlinearities, most notably the Cubic-Quintic (CQ) and Septic nonlinearities. I present, to the best of my knowledge, the first exact and stable solutions for these system. I then expand our method to include solutions for the Gross-Pitaevskii (GP) equation. In applying the F-expansion technique for this equation I encounter the Ricatti equation which cannot in general be solved for many systems. I present a large class of systems for which there are solutions to the Ricatti equation and from these solutions construct a large class of exact solutions for the GP equations. For constant strength of the external quadratic potential our solutions decay and we need an external gain to maintain these solutions, but for a sinusoidal form of the potential strength, one can obtain stable solutions. Additional solutions can be found for more complex forms of the diffraction coefficient and potential strength. xi I then study the effect of a linear potential on a NLSE with a cubic nonlinearity. The addition of the potential is studied in great detail and new solutions are obtained for this system. Finally, I analyze the stability of the aforementioned solutions. It is established that in most circumstances the obtained solutions are modulationally stable when dispersion management is used. Keywords Јacobi elliptic function nonlinear Schrödinger equation Gross-Pitaevskii Field Nonlinear optics Specialized Field Mathematical physics applied to nonlinear optics UDK number: 535.58 (043.3) Contents 1 Introduction 1 1.1 Motivation and aims . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Historical development of nonlinear optics . . . . . . . . . . . . . . . 2 1.3 Nature of the problem and theoretical methods employed . . . . . . . 4 1.3.1 Derivation of the Nonlinear Schro¨dinger Equation (NLSE) . 5 1.3.2 Terms of the NLSE . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Development of soliton solutions . . . . . . . . . . . . . . . . . . . . . 11 1.4.1 Vortex solitons . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.5 Methods of finding soliton solutions . . . . . . . . . . . . . . . . . . . 14 1.5.1 Inverse-Scattering theory . . . . . . . . . . . . . . . . . . . . . 15 1.5.2 Self-similar method . . . . . . . . . . . . . . . . . . . . . . . . 19 1.5.3 Hirota’s method . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.5.4 The Jacobi elliptic function expansion method . . . . . . . . . 21 1.6 Properties of Jacobi elliptic functions . . . . . . . . . . . . . . . . . . 24 1.7 Scope and structure of the Thesis . . . . . . . . . . . . . . . . . . . . 27 1.8 Published papers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2 F-expansion technique and its application to the general third- order NLSE 30 xii 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 F-expansion technique . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.3 Application of the F-expansion technique . . . . . . . . . . . . . . . . 36 2.3.1 Results for (3+1)-D NLSE with anomalous dispersion . . . . . 36 2.3.1.1 Case 1: Arbitrary β(z) and γ(z) . . . . . . . . . . . 38 2.3.1.2 Case 2: Arbitrary β(z) and χ(z) . . . . . . . . . . . 51 2.3.2 Results for (1+1)D NLSE . . . . . . . . . . . . . . . . . . . . 53 2.3.3 Results for (2+1)D NLSE . . . . . . . . . . . . . . . . . . . . 54 2.3.4 Results for (3+1)-D NLSE with normal dispersion . . . . . . . 55 2.4 Adaptation of the F-expansion technique for higher order nonlinearities 61 2.4.1 Results for n = 1 and n = 2 . . . . . . . . . . . . . . . . . . . 67 2.4.2 Results for n = 3 . . . . . . . . . . . . . . . . . . . . . . . . . 68 3 Application of the F-expansion technique to the Gross-Pitaevskii equation 70 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.2 Application of the F-expansion technique . . . . . . . . . . . . . . . . 71 3.3 Solutions for proportional α and β . . . . . . . . . . . . . . . . . . . 73 3.3.1 Solutions for constant α and β . . . . . . . . . . . . . . . . . . 75 3.3.1.1 Case 1: α and β of the same sign . . . . . . . . . . . 75 3.3.1.2 Case 2: α and β of opposite sign . . . . . . . . . . . 80 3.3.2 Solutions for sinusoidal α and β . . . . . . . . . . . . . . . . . 84 3.3.2.1 Case 1: α0 and β0 of the same sign . . . . . . . . . . 84 3.3.2.2 Case 2: α0 and β0 of the opposite sign . . . . . . . . 89 3.4 Other systems with solvable Ricatti equations . . . . . . . . . . . . . 90 xiii 3.4.1 Solution method . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.4.1.1 Case 1: A and B are constants . . . . . . . . . . . . 94 3.4.1.2 Case 2: A = 0 and B = B(x) . . . . . . . . . . . . . 95 3.4.2 Application of the solution method . . . . . . . . . . . . . . . 96 3.4.2.1 Example 1: β = 1 2 (e−δt + 1) . . . . . . . . . . . . . . 98 3.4.2.2 Example 2: β = ∑N n=0 βnt n . . . . . . . . . . . . . . 99 3.4.2.3 Example 3: β = β˜ ( 1− D B1t−B0 ) . . . . . . . . . . . . 102 4 Application of the F-expansion technique to the NLSE in a linear electric field 104 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.3.1 Case 1: Constant ǫ and β . . . . . . . . . . . . . . . . . . . . 107 4.3.2 Case 2: Constant ǫ and β = β0 cosΩt . . . . . . . . . . . . . . 108 4.3.3 Case 3: Constant β and ǫ = ǫ0 cosΩt . . . . . . . . . . . . . . 111 4.3.4 Case 4: β = β0 cosΩt and ǫ = ǫ0 cosΩt . . . . . . . . . . . . . 115 5 Two-component systems 117 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.2 Nonlinear Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.2.1 Application of the F-expansion technique . . . . . . . . . . . . 118 5.2.2 Case 1: f1 = ±g2 and f2 = ±g1. . . . . . . . . . . . . . . . . . 121 5.2.3 Case 2: c = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6 Stability Analysis 125 xiv 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.1.1 Modulational Instability . . . . . . . . . . . . . . . . . . . . . 127 6.2 Nonlinear Schro¨dinger equation and its transformation . . . . . . . . 128 6.3 Variational approach of the modulation stability . . . . . . . . . . . . 130 6.3.1 Case without chirp . . . . . . . . . . . . . . . . . . . . . . . . 132 6.3.2 Case with chirp . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.4 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 136 6.4.1 Split-step fast Fourier method . . . . . . . . . . . . . . . . . . 136 6.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 6.5 Analysis of stability for the Gross-Pitaevskii equation . . . . . . . . . 141 7 Conclusion 146 xv List of Figures 1.1 Jacobi elliptic functions as a function of u and M : (a) sn, (b) cn, (c) dn, (d) Jacobian amplitude am. . . . . . . . . . . . . . . . . . . . . . 26 2.1 Periodic traveling wave solutions as functions of the propagation dis- tance, for a0 = 0 (without chirp) and ǫ = 0. Intensity |u|2 for: (a) F = sn and (b) F = cn, presented as functions of k0x+ l0y+m0t and z. Coefficients: β(z) = cos(z), γ(z) = γ0 = 0, M = 0.9999, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0. . . . . . . . . . . . . . . . . . . . . 40 2.2 Periodic traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.1), except for a0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . 40 2.3 Combined intensity distributions of the periodic waves 1 and 4 from Table (1.1) (F = sn and F = ns), as functions of the propagation distance, with ǫ = 1 for: (a) a0 = 0 (no chirp) and (b) a0 = 0.1 (with chirp). Other parameters are the same as in Fig. (2.1). . . . . . . . . 41 2.4 Periodic traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.1)(a), except for: (a) M = 0.999, (b) M = 0.99, (c) M = 0.9 and (d) M = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.5 Solitary wave solutions without chirp. The setup and parameters are as in Fig. (2.1), except for M = 1. . . . . . . . . . . . . . . . . . . . . 42 2.6 Solitary wave solutions with chirp. Setup is the same as in Fig. (2.5), except for a0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 xvi 2.7 Combined intensity distributions of the solitary wave solitons 1 and 4 (F = sn and F = ns). The setup is as in Fig. (2.3) except forM = 1. 43 2.8 Solitary traveling wave solutions, as functions of the propagation dis- tance. The setup and parameters are the same as in Fig. (2.5), except for F = cn and: (a) β(t) = 1, b0 = 0.1, a0 = 0, (b) β(t) = 1, b0 = 0.1, a0 = 0.02 (c) β(t) = −1, b0 = 0.1, a0 = 0 and (d) β(t) = −1, b0 = 0.1, a0 = 0.02. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.9 The nonlinearity streingth given the integrability condition, as a func- tion of the the propagation distance z and the initial chirp a0. The parameters are the same as in Fig. (2.5), except for F = cn and: (a) N = 1 and (b) N = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.10 Solitary traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.5), except for: (a) F = sn, a0 = 0.2 (b) F = cn, a0 = 0.2, (c) F = sn, a0 = −0.1 and (d) F = cn, a0 = −0.1. . . . . . . . . . . . . . 46 2.11 Solitary traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.5), except for F = cn and: (a) b0 = 0 (b) b0 = 0, a0 = 0.1, (c) b0 = 0.5 and (d) b0 = −1. . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.12 Solitary wave solutions without chirp. The intensities |u|2 are plotted of solutions: (a) 1 (F = sn), (b) 2 (F = cn), (c) 7 (F = sc), (d) 8 (F = cs), (e) 4 (F = ns) and (f) 5 (F = nc). The parameters are otherwise the same as in Fig. (2.1), except M = 1 and γ(z) = −0.05. 48 2.13 Traveling wave soliton solutions without chirp. The parameters are the same as in Fig. (2.12), except M = 0.9999. . . . . . . . . . . . . 49 2.14 Solitary and travelling wave soliton solutions with chirp. The param- eters in (a) and (b) are the same as in Fig. 2.12 (a) and (b), and the parameters in (c) and (d) are the same as in 2.13 (a) and (b), except a0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 xvii 2.15 Phase of the solutions B as functions of the propagation distance and one transverse variable, assuming y = t = 0. The parameters are the same as in Fig. (2.5), except for F = cn and: (a) a0 = 0, and (b) a0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.16 Solitary wave solutions for β(z) = 1 + sin(z), f0 = 1, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. The remaining parameters are: (a) a0 = 0, F = sn and χ(z) = 1, (b) a0 = 0, F = cn and χ(z) = −1, (c) a0 = 0.02, F = sn and χ(z) = 1, and (d) a0 = 0.02, F = cn and χ(z) = −1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.17 Solitary wave solutions for β(z) = 1, χ(z) = −(2+sin(z)) and F = cn. Other parameters are the same as in Fig. (2.16), except for: (a) a0 = 0 and (b) a0 = 0.02. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 2.18 Values of γ(z) needed to achieve solution, as given in Eq. (2.59). Parameters in (a) are the same as given in Fig. (2.16), whereas the parameters in (b) are the same as given in Fig. (2.17). In both graphs a0 = 0 for the lower plot and a0 = 1 for the upper plot. . . . . . . . . 54 2.19 Solutions for a sinusoidal chirp of the (2+1)-D NSLE. The parameters are: α(z) = sin(z), γ(z) = sin(z), f0 = 1, b0 = 1, e0 = 0, k0 = l0 = 1, ω0 = 0, ǫ = 0 and: (a) M = 1, F = sn, (b) M = 1, F = cn, (c) M = 0.999, F = sn and (d) M = 0.999, F = cn. . . . . . . . . . . . . 56 2.20 Solutions for β(z) = − 1 2(1+z)2 , γ(z) = cos(z)/2, F = sn, M = 1. Other parameters are the same as in Fig. (2.19) except for: (a) α = 0 and (b) α = 1 + z. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.21 Solitary wave solutions as functions of the propagation distance, for a0 = 0 (without chirp) and ǫ = 0. Intensity |u|2 of: (a) solution 1 and (b) solution 2 from Table (1.1), presented as functions of k0x+ l0y+ m0t and z. Other parameters are: β(z) = cos(z), γ(z) = γ0 = −0.05, M = 1, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0. . . . . . . . . . . . . 59 xviii 2.22 Solitary wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.21), except for a0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.23 Traveling wave solutions without chirp. The setup and parameters are as in Fig. (2.21), except for M = 0.9999. . . . . . . . . . . . . . . 60 2.24 Traveling wave solutions with chirp. Setup is the same as in Fig. (2.23), except for a0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . 61 2.25 Soliton solutions for the cubic-quintic model n = 2 as a function of the propagation distance, for: (a) a0 = 0 (without chirp) and (b) a0 = 0.1 (with chirp) for the F = sech solution. Intensity |u|2 presented as a function of k0x + l0y − Ω0t and z. Other parameters are: β(z) = cos(z), γ(z) = γ0 = −0.05, M = 1, b0 = 1, e0 = 0, ǫ = 0, k0 = l0 = −Ω0 = 1, φ0 = 0. . . . . . . . . . . . . . . . . . . . . . . . . 67 2.26 Soliton solutions for the septic model as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.25), except for n = 3, F = tanh and ǫ = 1. . . . . . . . . . . . . . . . . . . 68 2.27 Periodic traveling wave solutions for the septic model as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.26), except for M = 0.99 and F = sn. . . . . . . . . . . . . 69 3.1 Decaying bent solitary wave solutions to GPE as functions of time, for b0 = 1. Intensity |u|2 for: (a) F = tanh and (b) F = sech presented as functions of k0x+ l0y +m0z and t. Other parameters are: β = 1, α = 1, γ(t) = −0.05, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. 76 3.2 Decaying straight soliton solutions to GPE as functions of time. The setup and parameters are the same as in Fig. (3.1) except for b0 = 0. 77 3.3 Decaying traveling wave solutions, given in terms of JEFs for: (a) F = sn and (b) F = cn. The setup and parameters are the same as in Fig. (3.2), except for M = 0.99. . . . . . . . . . . . . . . . . . . . 77 xix 3.4 Bent soliton solutions as functions of time for: (a) F = sn and (b) F = cn. The setup and parameters are the same as in Fig. (3.1), except for γ(t) = 3/ √ 2, the critical value of γ. . . . . . . . . . . . . . 78 3.5 Straight soliton solutions as functions of time. The parameters are the same as in Fig. (3.4) except for b0 = 0. . . . . . . . . . . . . . . . 79 3.6 Traveling wave solutions in terms of JEFs. The parameters are the same as in Fig. (3.5), except for M = 0.99. (a) F = sn and (b) F = cn. 79 3.7 Solitary wave solutions to GPE as functions of time. Values of pa- rameters are: β = 1, α = −1, γ(t) = −0.05, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. Intensity |u|2 for (a) F = tanh, b0 = 2 and for the remaining graphs F = sech with: (b) b0 = 2, (c) b0 = 0 and (d) b0 = −2 presented as functions of k0x+ l0y+m0z and t, where t is shown from 0 to T/4− 0.01. . . . . . . . . . . . . . . . . 81 3.8 Solitary wave solutions to the GPE as functions of k0x+l0y+m0z and t, where t is shown from 0 to T/2− 0.01. The setup and parameters are the same as in Fig. (3.7)(d) except for: (a) a0 = 0.1, (b) a0 = 0.5, (c) a0 = 1 and (d) a0 = 2. . . . . . . . . . . . . . . . . . . . . . . . . 82 3.9 Traveling wave solutions to the GPE as functions of k0x+ l0y +m0z and t, where t is shown from 0 to T/2− 0.01. The setup and param- eters are the same as in Fig. (3.7) except for M = 0.99, a0 = 1 and: (a) F = sn, b0 = 0, (b) F = cn, b0 = 0, (c) F = sn, b0 = 10 and (d) F = cn, b0 = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.10 Soliton solutions to Gross-Pitaevskii equation as functions of time, for the sine case: α(t) = α0 sin(Ωt), β(t) = β0 sin(Ωt). Intensity |u|2 presented as functions of k0x + l0y + m0z and t for: (a) F = tanh, b0 = 0, (b) F = sech, b0 = 0, (c) F = tanh, b0 = 1 and (d) F = sech b0 = 1. Other parameters are: α0 = 1, β0 = 1, Ω = 1, γ(t) = −0.05, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. . . . . . . . . . . . . 86 xx 3.11 Soliton solutions to Gross-Pitaevskii equation as functions of time for the cosine case: α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt). Other parameters are the same as in Fig. (3.10) . . . . . . . . . . . . . . . . 87 3.12 Traveling wave solutions to Gross-Pitaevskii solutions in terms of JEFs for the sine and the cosine case. The parameters for (a) and (b) are the same as in Fig. (3.10) (a) and (b), and the parameters in (c) and (d) are the same as in Fig. (3.11) (a) and (b), except for M = 0.99. For figures (a) and (c) we have F = sn and for figures (b) and (d) we have F = cn. . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.13 Soliton solutions to Gross-Pitaevskii equation as functions of time for the cosine case: α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt). Other parameters for figures (a) and (b) are the same as in figures Fig. (3.10)(c) and (d), respectively, except for b′0 = −b0 = 1. . . . . . . . 90 3.14 Graphs of parameters for Example 1: (a) α(t), (b) a(t) for a0 = 0 and (c) a(t) for a0 = 1, for δ = 0.01, 0.1, 1, 10 (top to bottom). . . . . 99 3.15 Intensity distribution |u|2 for the solution of Example 1. (a): no gain/loss. (b): γ = −0.05. Here F = sech. Other parameters are: a0 = f0 = k0 = l0 = m0 = ω0 = 1, b0 = ǫ = 0 and δ = 5. . . . . . . . . 99 3.16 Graphs of parameters for Example 2: (a) α(t), (b) a(t) for a0 = 0, and (c) a(t) for a0 = 1. Parameters: N = 0, 1, 2, 3, 4 (top to bottom at t = 0.5 for α, bottom to top at t = 3 for a), βn = 1. . . . . . . . . 100 3.17 Graphs of parameters for β(t) = cos(Ωt): (a) α(t), (b) a(t) for a0 = 0 and (c) a(t) for a0 = 1. Here Ω = 6, 7, 8, 9, 10. Curves with higher peaks correspond to lower values of Ω. . . . . . . . . . . . . . . . . . 100 3.18 Intensity distribution in Example 2, with b0 = 0 (left), and b0 = 5 (right). Here F = sech. In both cases, Ω = 8. Other parameters: a0 = f0 = k0 = l0 = m0 = ω0 = 1, γ = ǫ = 0. . . . . . . . . . . . . . . 102 3.19 Intensity distribution in case 3, with no gain/loss (γ = 0). Here, F = sech. Other parameters D = 10, B0 = −1, a0 = B1 = f0 = ω0 = e0 = k0 = l0 = m0 = β˜ = 1, and b0 = 0. . . . . . . . . . . . . . . . . . 103 xxi 4.1 :Soliton solutions for β and ǫ constant as functions of time. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) is presented as a function of k0x + l0y +m0z and t. Other parameters are: M = 1, β = 1, b0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0.2, δ = 0. Figures (a) and (b) are with no chirp or gain: a0 = 0, γ(t) = 0, Figures (c) and (d) with chirp: a0 = 0.1, γ(t) = 0 and Figures (e) and (f) are with chirp and gain: a0 = 0.5, γ(t) = 3/(2t). . . . . . . . . 109 4.2 Traveling wave solutions for β and ǫ constant as functions of time. The parameters are the same as in Fig. (4.1) except for M = 0.995. . 110 4.3 Soliton and traveling wave solutions for β = β0 cosΩt and ǫ constant. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) presented as a function of k0x+ l0y+m0z and t. Coefficients: β0 = 1, γ(t) = 0, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0.1, Ω = 1, γ(t) = 0, δ = 0. Figures (a), (b) depict solitary waves with no chirp: a0 = 0, M = 1, Figs. (c), (d) solitary waves with chirp: a0 = 0.1, M = 1, and Figs. (e), (f) traveling waves with chirp: a0 = 0.1, M = 0.99999. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.4 Soliton and traveling wave solutions for β constant and ǫ = ǫ0 cosΩt. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) presented as a function of k0x + l0y +m0z and t. Other parameters are: β = 1, γ(t) = 0, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ0 = 1, Ω = 1, δ = 0. Figures (a), (b) depict solitary waves with no b0: M = 1, b0 = 0, Figures (c), (d) solitary waves with b0: M = 1, b0 = 1, and Figures (e), (f) traveling waves with b0: M = 0.99, b0 = 1. 113 4.5 Soliton and traveling wave solutions for β constant and ǫ = ǫ0 cosΩt with chirp and gain. The parameters are the same as in Fig. (4.4) except for a0 = 0.5, γ(t) = 3/(2t). . . . . . . . . . . . . . . . . . . . . 114 xxii 4.6 Soliton and traveling wave solutions for β = β0 cosΩt and ǫ = ǫ0 cosΩt as functions of time. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) presented as a function of k0x + l0y + m0z and t. Coefficients: β0 = 1, γ(t) = 0, a0 = 0, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ0 = 3, Ω = 1, δ = 0. Figures (a), (b) depict solitary waves with no chirp: a0 = 0, M = 1, Figures (c) and (d) depict solitary waves with chirp: a0 = 0.1, M = 1 and Figures (e) and (f) depict traveling waves with chirp: a0 = 0.1, M = 0.999 . . 116 5.1 Traveling wave solutions for F = dn andG = nd constant as functions of time. Intensities: (a) |u1|2 and (b) |u2|2 are presented as a function of k0x and z for β(z) = β0 cosΩz. Other parameters are: M = 0.99, b0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, Ω = 1, β0 = 1, γ = 0, ǫ = φ = 1 and δ = −1. . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.2 Traveling wave solutions as functions of time. The parameters are the same as in Fig. (5.1) except for a0 = 0.1. . . . . . . . . . . . . . . 124 5.3 Traveling wave solutions as functions of time. The parameters are the same as in Fig. (5.1) except for F = sn and G = dc. . . . . . . . 124 6.1 Numerical simulation of the light bullet from Fig. (2.5)(b). Initial data from Eq. (2.56) are propagated according to Eq. (2.30) for 90 diffraction lengths along the z axis. Only the dependence on t is shown. The initial profile is noted by open circles. The curves to the left present intensity profiles at the left turning point, the curves to the right the profiles at the right turning point. The curves at the center are snapshots of the profiles passing approximately through the point t = 0 (i.e., the frames closest to t = 0 are recorded). Three sets of 15 profiles are overlapped at different z points, to show that no instabilities develop. . . . . . . . . . . . . . . . . . . . . . . . . . . 126 xxiii 6.2 (a) Nonlinearity parameter d for solutions cn and sn. (b) The growth rate parameter γ for dark an bright solitons as a function of K for the case σση = 1. Modulational instability occurs for values ofK depicted on the respective graphs. The solid lines represent the theoretical calculation of K using Eq. (6.33), and the square and circle dots are values of γ measured using numerical simulations, in which the dark and bright solitons, respectively, were perturbed by a small wave of the given wave number K. . . . . . . . . . . . . . . . . . . . . . . . . 132 6.3 (a) Theoretical values of log |U/U0| based on Eqs. (6.25)-(6.26) for a0 = 0. (b) Theoretical values of |U/U0| based on Eqs. (6.25)-(6.26) for a0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4 Maximum amplitude of oscillations plotted against Ω and K. Param- eters used: β0 = 0, β1 = 1, d = 8/3, (a) a0 = 0, (b) a0 = 0.02, (c) a0 = 0.1 and (d) a0 = 0.3. The contour boundaries, starting from the color purple are: 1, 2, 5, 10, 20, 50. . . . . . . . . . . . . . . . . . . . 137 6.5 Development of modulational instability for the bright soliton for three different values of z. Here, x is the direction of perturbation, y is the direction of the soliton and t is the remaining transverse direc- tion. Lower-wavelength colors (i.e. towards the color red) indicate a higher value of |u|2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.6 Development of modulational instability for the dark soliton for three different values of z. Here, x is the direction of perturbation, y is the direction of the soliton and t is the remaining transverse direction. Lower-wavelength colors (i.e. towards the color red) indicate a higher value of |u|2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.7 Evolution of the modulus of the perturbation amplitude |U | of bright solitons in time for different values of Ω: a) Ω = 1, b) Ω = 8. In both figures a0 = 0.3. Other parameters: K = 1, α0 = 0.3 [black, (upper) solid line]; K = 1, α0 = 0.1 [black, (upper) dashed line]; K = 4, α0 = 0.3 [red, (lower) solid line] and K = 4, α0 = 0.1 [red, (lower) dashed line]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 xxiv 6.8 Maximum amplitude of oscillations plotted against Ω and K. Pa- rameters used: β0 = 0, β1 = 1, d = 8/3, (a) α0 = 0.5, a0 = 0, (b) α0 = 0.05, a0 = 0 and (c) α0 = 0.05, a0 = 0.15. The contour boundaries, starting from the color purple, are: 2, 5, 10, 20, 50, 100. . 145 xxv List of Tables 1.1 Jacobi elliptic functions . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.1 Values of r-parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.1 Possible choices of F and G . . . . . . . . . . . . . . . . . . . . . . . 122 6.1 Stability cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 xxvi Chapter 1 Introduction 1.1 Motivation and aims The field of nonlinear optics (NLO) has made very important advances in the last 50 years and has yielded important theoretical and practical contributions to a di- verse range of fields: medicine, chemistry, engineering, telecommunications, biology, optics, as well as other branches of physics (1; 2). The biggest challenge in nonlinear optics is finding solutions to the nonlinear partial differential equations that naturally occur in the systems studied (3; 4). In particular, localized solutions to nonlinear equations that exhibit no change of form and which remain unperturbed are called solitons and have important applications in telecommunications, more specifically in the utilization of optic fibers. The field of finding exact solutions to nonlinear patrial differential equations (PDEs) is huge and many sophisticated methods have been devised to find them. Among the most prominent are: the inverse scattering technique (5; 6), the self- similar method (7), Hirota’s method (8), the Wronskian technique (9), the Back- lund transformation (10), the truncated Painleve´ expansion (11) and the function expansion method (12). The advent of powerful symbolic computation systems, such as Mathematica, has made new methods of finding direct solutions feasible. In particular, it became possible to effectively implement function expansions of the solutions in terms of 1 various functions, most notably the exponential functions, including the hyperbolic functions, and various elliptic integrals and functions, such as the Weierstrass elliptic functions and also the Jacobi elliptic functions (JEFs). The term F-expansion technique unless indicated otherwise will be used in this Thesis to refer to the expansion of solutions in terms of Jacobi elliptic functions. The functions themselves have properties which make them suitable for solving, with some modification, a large array of nonlinear partial differential equations. The main goal of this Thesis will be to expand the applicability of the F-expansion method to several important forms of the Nonlinear Schro¨dinger Equation (NLSE) and the Gross-Pitaevskii (GP) Equation. 1.2 Historical development of nonlinear optics Following (2), we will give a brief outline of the historical development of the field of nonlinear optics. The field of nonlinear optics began in the 19th century when John Kerr discov- ered the first instance of a quadratic electro-optic effect (13). Kerr observed the dependence of the index of refraction on the application of the electric field, a phe- nomenon that is now referred to as the DC Kerr effect. This is not to be confused with the optical Kerr effect, also known as the AC Kerr effect, where the change of the index of refraction comes purely from the light itself, though we shall for brevity sake refer to it as the Kerr effect in subsequent sections, as we will not be dealing with the DC Kerr effect in this Thesis. The discovery of the DC Kerr effect was immediately followed by the discovery of the linear electro-optic effect, i.e. the Pockels effect, named after its discoverer, Friedrich Pockels (14). Unfortunately, the field of nonlinear optics did not take off upon the discoveries of the Kerr and Pockels effects. Though the connection between light and electricity was well known, the theoretical basis for nonlinear effects in optics were still a mystery. However, with the advent of modern atomic theory in the early 20th century, i.e. the discovery of the true structure of the atom, it became possible through the effect of atomic polarization to develop he theoretical basis that would 2 explain these nonlinear effects. The first theoretical basis for a nonlinear optics phenomenon that did not in- volve an external electric field, the two-photon absorption, was developed by Maria Go¨ppert-Mayer in her Ph.D. Thesis (15). Go¨ppert-Mayer discovered that it is the- oretically possible for an atom to simultaneously absorb two photons leading to an energy excitation equal to the combined energies of the two photons. Like the Kerr effect, this effect depends on the intensity of light. Subsequently, Buckingham for- mulated the theory of molecular orientation in a strong magnetic field leading to a description of the optical Kerr effect (16) and then Piekara and Kielich developed the first theoretical formulation of the principles of nonlinear optics (17). Due to the high intensities of light required to observe the nonlinear effects, most nonlinear optics effects remain unobserved. This changed with the advent of the laser (18), allowing scientists to produce high-intensity coherent beams of light, after which the field of nonlinear optics underwent a period of rapid development. In particular, the field of experimental nonlinear optics was established with the first nonlinear effect observed, namely the second harmonic generation (SHG) (19). Var- ious other nonlinear optical effects were also very quickly observed: third harmonic generation (THG) (20), sum and difference frequency generation (SFG and DFG) (21), optical rectification (OR) (22), two-photon absorption (TPA) (23; 24), light- induced rotation (25; 26), stimulated Raman scattering (SRS) (27), the AC Kerr effect (28), self-phase modulation (SPM) (29; 30), cross-phase modulation (XPM) (31) and four-wave mixing (FWM) (32; 33). The discovery of these nonlinear effects has led to numerous practical applica- tion. One of the earliest was the so-called frequency doubling, as a result of SHG. This has allowed the construction of lasers emitting light in the visible range whose natural frequency, i.e. the difference in energy levels, would otherwise be somewhere in the infrared range (34; 35). Additionally, the effect of optical phase conjugation is used in holographic technology (36) and data encryption (37). NLO also has the prospect of being useful in the development of quantum computers, since the inter- action of two beams of light in a nonlinear medium has the possibility of producing a control-NOT gate, vital for the creation of a quantum computer (38). 3 One of the key applications of NLO is in telecommunication (2). Relatively quickly after the discovery of the first nonlinear effects, it was realized that optical wave-guides could be used to carry signals over a long range of distances (39). Gradually, through the improvement in the quality of the optical fibers (40; 41), the loss of signal was gradually reduced to levels where the use of fiber-optical cables for long-distance communication became feasible (42). The challenges today in NLO involve finding new materials that exhibit unique optical properties. Recently, a whole new class of materials, called metamaterials, was developed (43; 44). These materials exhibit striking properties such as a negative index of refraction (45; 46). As will be elaborated later on, these materials will allow us to produce stable solitary wave solutions for various physical systems, most notably the Nonlinear Schro¨dinger Equation (47). 1.3 Nature of the problem and theoretical meth- ods employed The Nonlinear Schro¨dinger Equation (NLSE) is one of the most prominent equations in nonlinear optics (48). It naturally arises in the study of the Kerr effect. Finding and classifying exact solutions to the equation remains an ongoing and active area of research within the field of NLO (49; 50). In systems governed by the NLSE, the index of refraction depends on the intensity of light this produces self-focusing, which, if balanced by the effects of diffraction or dispersion, can produce spatial or temporal soliton solutions to the NLSE, respectively (4). The solutions of the NLSE we obtain in this Thesis can most accurately be classified as spatio-temporal travelling wave and solitary wave solutions, though we will often refer to the solitary solutions as solitons. In this section, we will attempt to give a brief overview of the derivation of the equations which occur in this Thesis. All of the equations we will be dealing with ultimately stem from the NLSE. The derivation of the NLSE is a standard derivation in the field of nonlinear optics and here we will present a brief summary of this derivation, following (3), (4) and (51). 4 1.3.1 Derivation of the Nonlinear Schro¨dinger Equation (NLSE) The standard Maxwell’s Equations in differential form in the absence of external charge and current (in SI units) are: ∇×−→E = −∂ −→ B ∂t , (1.1) ∇×−→H = ∂ −→ D ∂t , (1.2) ∇ · −→D = 0, (1.3) ∇ · −→B = 0, (1.4) where −→ E , −→ H , −→ D , −→ B are, respectively, the electric field, the magnetic field, the electric displacement and the magnetic induction and ∇ = ∂ ∂x iˆ + ∂ ∂y jˆ + ∂ ∂z kˆ is the gradient operator. The relationship between the magnetic field and the magnetic induction is: −→ B = µ0( −→ H + −→ M), (1.5) while relationship between the electric field and the electric displacement is: −→ D = ǫ0 −→ E + −→ P , (1.6) where −→ M is the magnetic moment, which we assume to be equal to 0, and −→ P is the polarization response of the material. Here, ǫ0 and µ0 are coefficients of electric permittivity and magnetic permeability of vacuum, respectively. The polarization response −→ P can be given in the form of the following tensor: Pi(r, t) = ǫ0 ∞∑ k=1 χ (k) i,j1,...,jk k∏ l=1 Ejl , (1.7) where i, j1, . . . , jk are indices representing one of the three spatial directions, x, y or z, and Pi is the polarization response in the i direction, which we will separate into the linear −→ P (L) (k = 1) and the nonlinear −→ P (NL) (k > 1) responses. Assuming only the diagonal terms χxx to be present and mutually equal for all directions due to symmetry, we obtain −→ P (L) = ǫ0χxx −→ E . Plugging Eqs. (1.5)-(1.6) into Eqs. (1.1)- (1.4), then combining Equations (1.1)-(1.4) by applying the curl operator to Eq. (1.1) and plugging Eq. (1.2) into the RHS of Eq. (1.1), we obtain the standard 5 wave differential equation: ∇2−→E − n 2 c2 ∂2 ∂t2 −→ E = µ0 ∂2 ∂t2 −→ P (NL). (1.8) Here we have used the well-known identity: ∇×∇×−→E = ∇ · (∇ · −→E )−∇2−→E (1.9) and assumed ∇ · −→E = 0, which is only approximately true in the nonlinear regime. In practice, this approximation amounts to ignoring the longitudinal components of −→ E and −→ P (3). The parameter n given in the equation is the index of refraction defined by the formula n2 = 1 + χxx. If all terms in −→ P vanish except χ(1), then −→ E and −→ D are proportional, −→ P (NL) = 0 and one obtains all the results of standard linear optics. If, however, higher order terms exist, then one enters the field of so-called nonlinear optics. The primary difficulties that these higher order terms introduce is that they make the wave equations that are derived from Maxwell’s equations nonlinear. It can be shown, however, that in the presence of inversion symmetry many terms vanish, most notably the second-order terms (3). Of the third order terms of great interest is the term χxxxx which produces the so-called Kerr nonlinearity, which, as a result, makes the effective index of refraction dependant on the intensity of light. We now assume that the only term that exists is χxxxx and we assume that this term is real, i.e. we assume there is no nonlinear gain in our system. Switching to the frequency domain, we also assume the solution of the equation to be a wave envelope of the following form: −→ E = 1 2 (−→ A (x, y, z, t)ei(k0z−ω0t) + c.c. ) , (1.10) where ”c.c.” stands for the complex conjugate of the previous expression, ω0 is the central frequency and k0 = k(ω0) is the corresponding wave-number. The dispersion relation between parameters k and ω is k = ωn/c. The polarization term has the similar form, i.e: −→ P = 1 2 (−→p (x, y, z, t)ei(k0z−ω0t) + c.c.) . (1.11) 6 Plugging Eq. (1.10) into the third order term of Eq. (1.7) and neglecting third- harmonic generation we obtain −→p (NL) ≈ 3 4 ǫ0χxxxx|A|2−→A. (1.12) Applying the Fourier transform to Eq. (1.8) and taking into account the dispersion relation we obtain the following equation: (∇2 + k(ω)2)−→E (ω) = −µ0ω2−→P (NL)(ω). (1.13) Plugging in Eqs. (1.10)-(1.11) into Eq. (1.13) we obtain:( ∂2 ∂x2 + ∂2 ∂y2 + ∂2 ∂z2 + 2ik0 ∂ ∂z − k20 + k2(ω) )−→ A (x, y|ω) = −µ0ω2−→P (NL)(ω), (1.14) where −→ A (x, y|ω) is the Fourier transform of −→A (x, y, t). We now make the assumption that k is a slowly varying function of ω around ω0 and hence that we can write: k(ω) = ∞∑ j=0 kj j! (ω − ω0)j, (1.15) where k0 = n0ω0/c, k1 = dk dω |ω0 = 1vg is the inverse group velocity and k2 = d 2k dω2 |ω0 is the second order coefficient known as the group velocity-dispersion parameter (4). This term will be crucial in counteracting the Kerr effect and producing stable solitons. We will ignore terms higher than this order. Plugging in the truncated expression of Eq. (1.15) into Eq. (1.14) and applying the inverse Fourier transform, we obtain( ∂2 ∂x2 + ∂2 ∂y2 + ∂2 ∂z2 + 2ik0 ∂ ∂z − k20 + ( k0 + ik1 ∂ ∂t − k2 2 ∂2 ∂t2 )2)−→ A (x, y, t) = −µ0 ( ω0 + i ∂ ∂t )2−→p (NL). (1.16) We now apply the paraxial approximation which states that the envelope A of our beam is slowly varying in the longitudinal direction, in our case the z direc- tion. In other words, the following assumption hold on the derivatives of the wave amplitude A: ∣∣∣∣∣d 2−→A dz2 ∣∣∣∣∣≪ 2ik ∣∣∣∣∣d −→ A dz ∣∣∣∣∣ , ∣∣∣∣∣d 2−→A dx2 ∣∣∣∣∣ , ∣∣∣∣∣d 2−→A dy2 ∣∣∣∣∣ . (1.17) We also make a similar assumption for time, i.e. that the amplitude varies slowly in time. In practice, this will mean neglecting time derivatives higher than 2 on the 7 LHS and 1 on the RHS. We will also assume that −→ A maintains polarization and, as mentioned before, that its longitudinal component is negligible, in which case we can switch to the scalar form of our equations and drop the vector signs. Applying all these approximations, and plugging in Eq. (1.12), we obtain: ∂2A ∂x2 + ∂2A ∂y2 + ∂2A ∂z2 +2ik0 ∂A ∂z +2ik0k1 ∂A ∂t − (k21 + k0k2) ∂2A ∂t2 = −3 4 µ0ω 2 0ǫ0χxxxx|A|2A. (1.18) We now introduce the following change of variables: τ = t−k1z. This is time in the frame moving with the group velocity of the pulse. The slowly-varying-in-time approximation for the amplitude means that: k1 ∣∣∣∣∂A∂τ ∣∣∣∣≪ k0 (1.19) and, in accordance with the chain rule, the derivatives in the τ and z directions now become: ∂ ∂z → ∂ ∂z − k1 ∂ ∂τ , (1.20) ∂2 ∂z2 → ∂ 2 ∂z2 − 2k1 ∂ ∂z ∂ ∂τ + k21 ∂ ∂τ , (1.21) ∂ ∂t → ∂ ∂τ . (1.22) Combining Eq. (1.17) and Eqs. (1.19)-(1.29) into Eq. (1.18) and dividing the expression by 2k0, we obtain i ∂A ∂z + 1 2k0 ( ∂2A ∂x2 + ∂2A ∂y2 ) − k2 2 ∂2A ∂τ 2 = −3µ0ω 2 0ǫ0χxxxx 8k0 |A|2A. (1.23) We can now perform a re-scaling of τ to align it with the other transverse variables: τ ′ = τ( √|k0k2|). From now on, we will for simplicity write τ ′ as t (we will not be using the original “real” time anymore in our Thesis). We will also re-scale the amplitude of the beam so that the integral squared is equal to the power of the beam, rather than the intensity, by introducing a new variable: u = A √ Seff , where Seff is the effective area of the beam. After re-scaling Eq. (1.23), we finally obtain: i ∂u ∂z + β(z) 2 ( ∂2u ∂x2 + ∂2u ∂y2 + s ∂2u ∂t2 ) + χ(z)|u|2u = iγ(z)u, . (1.24) This equation is known as the the time-dependent Nonlinear Schro¨dinger Equation (NLSE) with cubic nonlinearity (the order of u in the nonlinear term). Unless 8 otherwise indicated, we will assume all our NLSEs to have cubic nonlinearity. The form of the NLSE given in Eq. (1.24) will serve as the basis of the remainder of this Thesis. Here, β(z) = 1 k0 (1.25) is the diffraction/dispersion parameter, χ(z) = 3µ0ω 2 0ǫ0χxxxx 8k0Seff (1.26) is the strength of nonlinearity and s = −sgn(k0k2) is the sign parameter which describes whether the dispersion is normal s = −1 or anomalous s = 1. We have also introduced a new parameter γ(z) which will describe the external gain (for positive γ) or loss (negative γ) given to the system. Another thing to note about Eq. (1.24) is that we have purposely chosen all our parameters to be dependent on the longitudinal variable z. Furthermore, we allow all our parameters, especially β, to be both positive and negative. It is possible for β to be negative when the index of refraction is negative, which has been achieved for various types of meta- materials (44). Our goal is to use the variability of these parameters, especially the regimes where β oscillates between positive and negative values, to obtain novel spatio-temporal travelling wave solutions to the multidimensional NLSE. 1.3.2 Terms of the NLSE In this subsection we examine the appearance of various nonlinear terms in the NLSE which will be used throughout this Thesis. An important fact about Eq. (1.24) is that it is extremely versatile and can be adapted for a whole range of systems. We have already seen that the equation covers both the normal and anomalous dispersion for s = −1 and s = 1, respectively. However, we can also obtain the time-independent NLSE that ignores the effect of temporal dispersion if we put s = 0. In addition, we can place our beam within waveguides that support only a single mode, i.e. the fundamental mode, depending only on a given transverse direction. In this case, we can factorize the dependence of u along that direction out of our NLSE, and in this case Eq. (1.24) no longer contains the selected transverse derivative (4). We can thus write the most general 9 form of the NLSE as follows: i∂zu+ β(z) 2 N∑ i=0 si∂ 2 xi u+ χ(z)|u|2u = iγ(z)u, (1.27) where x0 is time and xi, i > 0, are the spatial transverse variables. Here s0 = −1, 0, 1, depending on whether we are looking at the case of normal dispersion, the time-independent equation, or anomalous dispersion, respectively, and si = 0, 1 for i > 0 depending on whether there is a waveguide perpendicular to the xi direction or not, respectively. Obviously, N = 1, 2. Note that in Eq. (1.27) the partial derivatives are written as indices, a notation we will frequently use. We will also generalize the NLSE to higher-order nonlinearities. Instead of looking at only the third-order nonlinearity, we will make the following substitution in Eq. (1.24): χ |u|2 u→ χ1 |u|2 u+ χ2 |u|4 u+ · · ·+ χn |u|2n u. (1.28) We now turn to the other terms that can be added to our equation. The first term that we will add to Eq. (1.24) is the external potential V (x, y, z, t)u. The potential is typically added for equations in which the time t is the longitudinal variable and all the spatial variables are transverse. In other words, t and z exchange roles. In that case, it is convenient to study potentials of the form V (x, y, z, t) = f(t)(xp+yp+zp), where f is some function of time. For p = 2 one obtains the Gross-Pitaevskii equation in the standard harmonic potential, and for p = 1, one obtains the NLSE in a linear potential. Both of these cases will be covered in this Thesis. Finally, we will also examine the additional optical effects. In our derivation of Eq. (1.24) we assumed only the term χxxxx exists, which leads to the optical effect known as self-phase modulation (SPM) that is accounted for in Eq. (1.24). However, other terms may also exist if the beam is not linearly polarized or is indeed itself a combination of two co- or counter-propagating beams. In this case, other components of the χ(3) can then cause the two beam components to interact with each other and, as a result, we can in many cases observe the effects of cross-phase modulation (XPM) and four-wave mixing (FWM). Accounting for these effects in a more general calculation than one done in Sec. 1.3.1, we obtain an equation in which the following substitution is done with respect to Eq. (1.24): ∣∣u2∣∣ u→ ∣∣u2∣∣u+ c ∣∣v2∣∣ u+ fv2u∗, (1.29) 10 where: c = 2 ( χ (3) xxyy(ω;ω, ω,−ω) + χ(3)xyyx(ω;ω,−ω, ω) + χ(3)xyxy(ω;ω, ω,−ω) ) 3χ (3) xxxx(ω;ω, ω,−ω) ,(1.30) f = χ (3) xxyy(ω;−ω, ω, ω) + χ(3)xyyx(ω;ω, ω,−ω) + χ(3)xyxy(ω;ω,−ω, ω) 3χ (3) xxxx(ω;ω, ω,−ω) , (1.31) v is the second component and u∗ is the complex conjugate of u (3). Here we have used the standard notation in NLO for the frequency, (ω;ω1, . . . , ωk), signifying that the nonlinearity coefficient χ corresponds to a wave of frequency ω = ω1 + · · ·+ ωk generated by waves with frequencies ω1, . . . , ωk, where negative frequencies correspond to the complex conjugate terms. The coefficients c and f correspond respectively to XPM and FWM. A symmetric equation can also be written for v, thus giving us a system of two equations with two unknowns. In our Thesis we will neglect the effects of FWM, i.e. we will have f = 0 in all our systems. This occurs when the mechanism through which nonlinearity occurs is slow, e.g. through photorefractive effect or the change of the temperature of the material (3). On the other hand, the CPM coefficient c can have a range of values, both arbitrarily large (52) and even negative (3; 53). For c = 1 we have the so-called Manakov system, which occurs for two mutually incoherent beams (54), whereas c = 2 occurs for coherent polarized beams in highly birefringent media (54; 55). In Sec. 5 we will obtain solutions for the two-component system when c = 3. 1.4 Development of soliton solutions In this section we will give a brief historical overview of the development of soli- tary wave, i.e. soliton, solutions to nonlinear partial differential equations, as well as methods developed for finding such solutions. It is typical in most literature on nonlinear optics to refer to all solitary wave solutions as solitons (4) and, as mentioned before, we will employ this terminology in this Thesis. Solitary waves were originally observed by John Scott Russell on the Union Canal at Hermiston in Edinburgh (56). Observing solitary waves travelling along the narrow canal Russell postulated that these solitary waves display many properties 11 of a classical particle. It is only in the middle of the 20th century that his discovery began to be developed. In the paper by Fermi, Pasta and Ulam (57) evidence first emerged that vibrating modes on an anharmonic lattice exhibit quasi-periodic behavior. Further research by Zabusky and Kruskal confirmed this (58). In their study of a variant of the Korteweg-de-Vries (KdV) equation, it became clear that several distinct waves form in the waves which seem to retain their identity even after colliding and interacting with the other waves. This led to the coining of the term “soliton” emphasizing the particle-like nature of these waves. Other systems in which such behavior was analyzed were the NLSE and the sine-Gordon (SG) equation. Solitons can be classified based on their shape as either bright solitons, where a signal of light travels on a dark background, dark solitons where a dark signal (absence of light) travels against a background of light, and so-called gray solitons, which are similar to dark solitons except that the intensity of light never drops to 0 as is the case of dark solitons (4). The first optical pulses in fibres, i.e. temporal solitons, for systems with Kerr nonlinearity were theoretically constructed by Hasegawa and Tappert and then ver- ified using computer simulations (59; 60). They proved that the dispersion term which causes pulses to widen could be balanced by the nonlinear Kerr effect which causes pulses to contract, giving us stable soliton solutions to the NLSE (61). This was experimentally verified by Mollenauer (62) and the transmission of dark solitons was verified by Emplit and his group (63). Since then, optical pulses have been ex- tensively used in telecommunication for transferring clear signals that do not decay in time and also do not mix and interfere with each other (42; 64). The first evidence of spatial solitons in nonlinear optics was demonstrated by Bjorkholm and Ashkin (65), who used a dye laser passing through a chamber filled with sodium vapor to achieve this. Soon, spatial solitons were observed in other kinds of media (66), most notably within a semi-conductor wave-guide (67). A novel medium in which spatial solitons have been achieved are nematic crystals (68) and the solitons produced there are known as nematicons (69). The most attractive feature of nematicons is that due to the reorentational nonlinearity of 12 nematic crystals they can be generated at a relatively low power (70). The study of the propagation of light inside nematic crystals remains an active area of research (71; 72). Of all forms of solitons, the most complex form are the so-called vector solitons (73). They can be viewed as solitons which are elliptically polarized. These types of solitons typically arise in systems governed by Manakov equations (54). The two components may be either dark or bright, and indeed it is even possible to form a bright-dark soliton (74). The first experimental evidence of vector solitons was given by Cundiff and his group (75), while so-called phase-locked vector solitons were discovered by Tang (76). 1.4.1 Vortex solitons A subclass of solitons of particular interest are the so-called vortex solitons (77). These solitons arise naturally in media with cylindrical symmetry. The solution to such systems is proposed to be of the form: u(r, θ) = U(r)eiφ(θ), (1.32) where r = √ x2 + y2 is the transverse radius, U the amplitude and θ the azimuthal phase. The quantity: S = 1 2π ∮ dφ (1.33) is called the topological charge of the vortex, and is necessarily an integer. Typically, we have φ = mθ, in which case S = m. For S = ±1 one can obtain stable vortex solitons, while those with |S| > 1 are typically unstable and break-up into components. Vortex solitons arise when there is a singularity in the phase at the origin, i.e when S 6= 0. This means that the amplitude U must also be equal to 0 at the center of the vortex. Fundamentally, there are two different types of vortex solitons: those arising in self-focusing media that support bright solitons, and those arising in self-defocusing media that support dark solitons. In the case of vortex solitons in self-focusing media one obtains a doughnut-shape vortex soliton. These are usually highly unstable and split up into components due to azimuthal instability (78), and these components 13 subsequently fly away from each other (79). Several proposals have been made for stabilizing vortex solitons in self-focusing media, most prominent being to introduce competing nonlinearities such as in the cubic-quintic model (80) and to introduce a non-local nonlinear response (81). Stable solitons were observed for the first time in (79) where partially incoherent light was used to dampen azimuthal instability and the theory the existence of this stability was developed in (82). In the case of dark solitons one obtains a topologically charged hole on a uniform background. Such solitons were first observed by Swartzlander and Law (83) and subsequently observed by Ref. (84) when dark soliton stripes undergo break-up due to transverse modulational instability. One can factor out the background from the self-defocusing NLSE with a suitable transformation (77) to obtain: i ∂u ∂z + 1 2 ∇2u+ (1− |u|2)u = 0. (1.34) Assuming the form of the solution to be as given in Eq. (1.32) with φ = mθ, one obtains the following amplitude equation: ∂2U ∂r2 + 1 r ∂U ∂r − m 2 r2 U + (1− U2)U = 0, (1.35) with boundary conditions U(0) = 0 and U(∞) = 1. Though no analytical solutions to this equation exist, it can be solved numerically (85). The study of vortex solitons remains a very active area of nonlinear optics that has branched out in several directions. In particular, modulation in the azimuthal direction began to be considered with the establishment of azimuthons (86; 87; 88). These structures proved very suitable for the study of vortices on lattices, both square (89) and hexagonal (90). Another approach has been the study of interaction of vortices with ordinary beams (91) and in general combining two beams to obtain rotating structures (92). Particularly useful has been the study of the interaction between two counter-propagating beams (93; 94; 95) where a range of stable rotating dipoles, tripoles, quadrupoles and other structures have been discovered (96). 1.5 Methods of finding soliton solutions A range of general techniques were developed for finding exact soliton solutions for large classes of partial nonlinear differential equations. Each technique has its own 14 scope and range of applicability. Also, some techniques do not give exact solutions, but merely approximations of solutions. Here we will only give a brief review of some of the most important techniques. 1.5.1 Inverse-Scattering theory The inverse scattering theory (IST) is a powerful technique whose primary goal is to obtain multi-soliton solutions for an arbitrary initial condition, using the formalism of quantum mechanics. The primary goal of the IST is to re-cast the given nonlinear (NL) equation into a linear Schro¨dinger equation (SE): ∂2xψ + (λ− u)ψ = 0, (1.36) where the solution to the NL equation u appears as the potential of the SE. The solutions therefore are to satisfy the Faddeev condition:∫ ∞ −∞ (1 + |u(x)|)|u(x)|dx <∞, (1.37) which is stronger than the normalization condition. Solving the SE, i.e. the so- called direct scattering problem, one obtains a number of parameters called the scattering data. Among these the most fundamental are the transmission T (t) and two reflection coefficients R(t) and L(t) (right and left) that arise out of considering the behavior at +∞ and −∞ of the continuous spectrum of eigenfunctions with positive eigenvalues, in particular of the so-called Jost solutions (97). As for the negative eigenvalues, one obtains a discrete spectrum of N eigenvalues, each with two normalization coefficients defines through the two Jost solutions. However, it turns out these two coefficients are related for each eigenvalue and we can thus only consider one set of coefficients cn. The IST needs to satisfy a range of conditions in order to work. The most important criterion is that of integrability. This means that the system must be chosen so that the eigenvalues of the SE remain constant (98). The criterion for determining integrability was devised by Lax (99). Lax formalism proposes a set of 15 equations satisfying: Lψ = λψ, (1.38) Mψ = ∂tψ, (1.39) where L and M are operators known as the Lax pair and indexes denote partial derivatives. It is then possible to combine the two expressions under the condition that ∂tλ = 0 to obtain the following equation: ∂tL+ [L,M] = 0, (1.40) where brackets denote the commutator of two operators. Thus, if the differential equation can be written in this form for some two operators L and M, then the eigenvalues of the system governed by L remain constant. In the standard IST L is taken to be the linear Schro¨dinger equation: L = −∂2x + u. (1.41) The disadvantage of this method is that there is no systematic way to identify whether a Lax pair exists for a differential equation, indeed it is even possible to find multiple pairs for a single equation. Another complicating feature is that the operators need not be scalar. They can be matrix operators of arbitrary rank which extends the usefulness of the IST to the Nonlinear Schro¨dinger equation. In the next step, one allows the scattering data to vary across time. The formulas that can be derived are usually straightforward provided that the eigenvalues do not change with time (98). Finally, one obtains the value of the potential to the SE u(x, t) at this new point in time, based on the scattering data, i.e. one solves the inverse scattering problem. In other words, the goal is to reconstruct the new form of the potential based on the new scattering data. This is achieved by solving the so-called Marchenko equation, also known as the Gel’fand-Levitan-Marchenko (GLM) equation. The GLM equation typically contains Ω, the so-called Marchenko kernel, obtained from the scattering data, and the unknown parameter of the equationK is directly related to u. It is worth noting that in general the GLM equation cannot be solved in closed form. However, when the form of F is separable, which occurs when the potential 16 is reflectionless, then solutions of arbitrary order can be found. This occurs when u(0, x) = Asech x. For A = N(N + 1) for some positive integer N one obtains an exact solution of N solitons, while for (N +1)(N +2) > A > N(N +1) one obtains in addition a travelling wave in the −x direction (98). In a paper by Zakharov and Shabat, the inverse scattering method was used to derive a general solution for the (1+1)D NLSE (100). This scheme was later generalized and refined by Ablowitz, Kaup, Newell and Segur (101) and is now known as the ZS-AKNS scheme. Zakharov and Manakov determined that the NLSE in (1+1)D is integrable (100). Instead of scalar operators, matrices were used along with two coupled differential equations. Zakharov and Shabat used the following system: Ψx = LΨ, (1.42) Ψt = MΨ, (1.43) where: Ψ =   Ψ1 Ψ2   , (1.44) L =   −iζ u −u iζ   , (1.45) M =   −2iζ2 + i|u|2 2ζu+ iux −2ζu+ iux −2iζ2 + i|u|2   , (1.46) the overline denotes the complex conjugate, ζ is called the spectral parameter, and the function u satisfies the compatibility condition uxt = utx only if u is a solution to the (1+1)D NLSE: iut + uxx + 2|u|2u = 0. (1.47) A similar set of formulas to Eqs. (1.45)-(1.46) for the equation iut−uxx−2|u|2u = 0 can be found in Ref. (98). The matrices L and M satisfy the so-called AKNS condition: Lt −Mx + [L,M] = 0, (1.48) which is similar to the condition given in Eq. (1.40). The standard method of discovering the form of M in the AKNS scheme is to assume that each entry in 17 M is an arbitrary polynomial in ζ and obtain a system of differential equations by plugging in this form into M (102). The system is underdetermined due to symmetry and therefore one may impose the condition of the NLSE to be solved for u (98). The key complication in the application of the ZS-AKNS scheme for the NLSE is that bound-state eigenvalues are complex (with a positive imaginary value) as opposed to purely imaginary (103). Indeed, the term ζ corresponds to √ λ in Eq. (1.36). Also, the bound-state eigenvalues can have a multiplicity larger than 1 thus we obtain nj coefficients cl,j , 0 ≤ l ≤ nj − 1 for each ζj ∈ C+ (103). We summarize here the most basic results obtained for the NLSE, following (103) and (97). The scattering coefficients acquire the following dependance for the NLSE: T (t) = T (0), (1.49) R(t) = R(0)e4iζ 2t, (1.50) L(t) = L(0)e−4iζ 2t, (1.51) cj(t) = cj(0)e −4iA2j t, (1.52) where: cj(t) = [ cnj−1,j(t) cnj−2,j(t) . . . c0,j(t) ] (1.53) and: Aj =   0 −1 0 0 · · · 0 0 −1 0 · · · 0 0 0 −1 . . . 0 0 0 0 . . . ... ... ... ... . . .   − iζjI (1.54) is an nj × nj-sized matrix (104). For nj = 1, one obtains for the lone normalization coefficient: c0,j(t) = cj(t) = cj(0)e 4iζ2j t. (1.55) The Marchenko equation for the NLSE has the following form: K(x, y, t)− Ω(x+ y, t) + ∫ ∞ x dz ∫ ∞ x dsK(x, s, t)Ω(s+ z, t)Ω(z + y, t) = 0, (1.56) 18 where the formula for the Marchenko kernel is given as: Ω(x, t) = 1 2π ∫ ∞ −∞ dζR(ζ, t)eiζx + N∑ j=1 nj−1∑ j=0 cl,j xl l! eiζnx. (1.57) Finally, after solving the Marchenko equation for K one obtains the solution u via the following formula: u(x, t) = −2K(x, x, t). (1.58) As was mentioned before, for R = 0 the problem becomes analytically solvable. The primary solutions obtained using the IST have been the N soliton solutions with the reflectionless potential as the initial condition by Ref. (6). These solutions are valid for a set of N distinct eigenvalues and are generalized to the case of eigenvalues of multiplicity larger than 1 in Ref. (97). 1.5.2 Self-similar method Self-similarity is an important phenomenon common in nature. It is also routinely studied in mathematics, most importantly in fractals, which occur in a wide range of physical contexts, including solitons (105). In the context of differential equations, the so-called self-similar method involves the reduction of variables of a nonlinear partial differential equation (7). This is achieved by utilizing a group of transformations that leaves the partial differential equation invariant (106). For example, for the 1D nonlinear Schro¨dinger given in Equation (1.47) the transform leaving the equation invariant is: u(x, t)→ ǫu(ǫx, ǫ2t). (1.59) Therefore, the self-similar form of the solution will be: u(x, t) = 1√ t U( x√ t ). (1.60) The variable X = x√ t , in this example, is called the self-similar variable. A well- known special case of self-similar solutions are the so-called travelling-wave solutions, where the similarity variable is taken to be of the form X = kx+ ωt for some non- zero k and ω. This can easily be generalized to a higher number of variables, which is what we will do in this Thesis. 19 One of the most practical applications of the self-similar method lies in trans- forming nonlinear partial differential equations (NLPDEs) with multiple transverse variables into NLPDEs with just a single transverse variable. In (107) and (108), self-similar solutions to the NLSE were found using Hermite polynomials. In (109), self-similar solutions are obtained for a NLSE with an arbitrary streingth of nonlin- earity and external potential quadratic with respect to the transverse variable. In (110), a generalization of the solutions found in (111) for the higher order nonlinear Schro¨dinger equation is provided using the self-similar method. An example of a chirped self-similar soliton to the generalized cubic-quintic NLSE is given in (112), although these solutions are only valid for the constant value of the chirp. 1.5.3 Hirota’s method Hirota’s method is an ingenious method for obtaining exact multi-soliton solutions to large classes of differential equations (8). The main idea of Hirota’s method is to reduce the differential equation to one or two dependent variables, traditionally labelled F and G, such that the differential equation becomes quadratic. For the NLSE, the proper choice is: u = G F , (1.61) where, without loss of generality, G is complex and F is real. The newly obtained equation also has to be expressible via the so-called bilinear form. The definition of the bilinear form or two functions is: DkxD l yD m z D n t F ·G = ( ∂ ∂x − ∂ ∂x′ )k ( ∂ ∂y − ∂ ∂y′ )l( ∂ ∂z − ∂ ∂z′ )m (1.62)( ∂ ∂t − ∂ ∂t′ )n F (x, y, z, t)G(x′, y′, z′, t′) |x=x′, y=y′, z=z′, t=t′ . The form resembles the product rule for derivatives, except that it is anti-symmetric. Occasionally, as is the case with the NLSE, this requirement splits the equation into a system of two or more equations, due to the fact that one has to group the anti- symmetric terms together with certain common factors that may differ for each group. For example, for the 1D-NLSE with constant coefficients, i ∂u ∂t + ∂2u ∂x2 + χ|u|2u = 0, (1.63) 20 one obtains plugging in Eq. (1.61) the following pair of equations: (iDt +D 2 x)G · F = 0, (1.64) D2xF · F = χ|G|2. (1.65) Care needs to be taken for the obtained system of equations not to be overdetermined (113). After obtaining a set of differential equations using Hirota’s bilinear form, one then takes the solutions to be of the form where F = ∑N i=0 fiǫ i, G = ∑N i=0 giǫ i and ǫ is some small parameter. For the NLSE, one can neglect the odd terms for G and even terms for F . Plugging in this form of the solution into the equation, expanding it, taking care to multiply out all the fractions, and collecting the terms according to the order of ǫ, we obtain a series of differential equations, one for each order of ǫ. The convenient form of an ansatz for functions fi and gi is the sum of exponentials. By determining where the series truncate, one can obtain soliton solutions of an arbitrary order. The parameter ǫ is then taken to be 1. Hirota’s method has been extensively applied to find multi-soliton solutions in various systems: the NLSE (114), the Gross-Pitaevskii equation (115) and others (116). 1.5.4 The Jacobi elliptic function expansion method The functional expansion (F-expansion) method is probably one of the simplest methods conceptually for solving nonlinear PDEs and the approach we will be using in this Thesis. The idea is to assume a solution of a form involving the sum of certain convenient set of functions called the expansion functions multiplied by a set of variables, i.e. coefficients, which have to be determined from the equation. Well- known examples include the polynomial expansion and the exponential function expansion. Once the form of the solution is plugged into the equation, the terms can be grouped by the expansion functions themselves. Since the expansion functions are non-zero, the sum of terms in each group has to be 0. This condition yields a set of algebraic or ordinary differential equations involving the unknown coefficients which can hopefully be solved. Frequently, due to the complexity of obtained systems, 21 Wu’s elimination method is applied (117; 118). Extensive algebraic calculations are usually performed with the aid of a symbolic computation package such as Maple or Mathematica. In this Thesis, we shall exclusively use Mathematica for such computations, and also for the plotting of all graphs. The drawbacks of the F-expansion method are obvious. The effectiveness of the method largely rests on the choice of expansion functions. A poor choice of expansion functions will result in a self-contradictory system, or one possessing only the trivial solution wherein all the coefficients are equal to 0. Also, even in the process of obtaining solutions one will seldom obtain the most general solutions, but instead will obtain one or more constraints that the coefficients of the nonlinear PDE must satisfy. Still, if the expansion functions and the general form of the solution are suitably chosen, one will obtain large classes of exact solutions that can subsequently be studied in experiments or in computer simulations. The Jacobi elliptic functions (JEFs) are a particularly convenient choice for the F-expansion technique. The JEFs are a generalization of the trigonometric functions that satisfy a second order NL differential equation in which the cube of the original function appears, a form particularly convenient for the standard third- order NLSE that we will be studying in this Thesis. In the appropriate limits the JEFs converge to both the elliptic and hyperbolic trigonometric functions allowing us to simultaneously find both travelling wave and solitary solutions. The JEF expansion method was first fully developed in Ref. (12). In this paper, the JEF expansion method was proposed as a generalization of the tanh-expansion method (119). It was applied to the modified KdV (mKdV) equation and the nonlinear Klein-Gordon equation and, as a result, exact solutions to these equations were obtained. A key concept developed in the paper was to use the form of the differential equation for the JEFs to match the highest order of the derivative terms with the term containing the nonlinearity. The JEF expansion method was further developed in (120) where a modified ansatz was used to find the solutions of the mKdV equation and a system of variant Boussinesq equations. Once developed, the JEF expansion method was extensively used for a vari- ety of systems. In particular, in a paper by Yan (121) an extensive theory is 22 developed for both travelling and non-travelling wave solutions and solutions are obtained for the general KdV-mKdV equation (encompassing both equations as special cases), the modified Kadomtsev-Petviashvili (mKP) equaton, the coupled Davey-Stewartson (DS) equation and the modified Boussinesq equation. The JEF expansion method has also been applied to the Benjamin-Bona-Mahoney (BBM) equation (122), the Witham-Boer-Kaup equation (123), the coupled KdV equations (124), Klein-Gordon-Schro¨dinger (KGS) equation (125), the higher-order nonlinear Schro¨dinger equation (126), the dispersive long wave equation (127; 128; 129), as well as many other systems (130; 131; 132; 133; 134; 135). The JEF expansion method has also been applied extensively to the Nonlinear Schro¨dinger equation. In (136), the NLSE emerges as a special case of a coupled system whose solutions were found using the JEF expansion technique. Solutions to the NLSE have also been given in Ref. (137) and Ref. (138) for a (1+1)-D NLSE with constant coefficients. In (139), a novel approach is taken in that the solution to the (2+1)-D NLSE is decomposed via its real and imaginary parts, rather than the amplitude and phase. The coefficients of the NLSE, except for the gain/loss coefficient, are also assumed to be constant in that paper. The solutions to the NLSE using the JEF expansion technique all had several drawbacks in common before the work done in this Thesis. All of them were usually solutions to the NLSE with constant coefficients instead of distributed coefficients. As mentioned before, with the invention of metamaterials there is considerable in- terest in nonlinear systems that are more general than the NLSE with constant coefficients. The second big disadvantage is that none of the obtained solutions ad- dress the problem of chirp. The phenomenon of chirp was given extensive treatment for the first time in (140) with regards to the NLSE and it is of considerable interest to include those results in the framework of the JEF expansion method. The final limitation concerned the dimensions of the solved NLSEs. A simple approach was needed that would produce solutions for the (2+1)-D and (3+1)-D NLSEs. The goal of this Thesis will be to provide a general ansatz for the JEF expansion method that surpasses these limitations and can be generalized to a wide range of systems. 23 1.6 Properties of Jacobi elliptic functions In this section a brief overview of Jacobi elliptic functions is given. The standard reference for the Jacobi elliptic function is Ref. (141). Jacobi elliptic functions (JEFs) can be viewed as generalization of the trigono- metric functions. The trigonometric functions are defined on a unit circle, defined by the equation x2 + y2 = 1, as the (x, y) coordinates of the point which is at an angle θ away from the point (1, 0), i.e. sin θ = y and cos θ = x. In the case of the Jacobi elliptic functions, however, we will use the ellipse whose equation is x2 + y 2 a2 = 1, where a > 1. The eccentricity of the ellipse is given as k = e = √ 1− 1 a2 . When a = 1, the eccentricity is 0 and in this case the ellipse becomes a unit circle. Given polar coordinates, the argument of the ellipse will be defined as u = ∫ θ 0 rdθ. It is important to note that the argument is not equal to the arc length in general, as it doesn’t take into account the change in radius of the ellipse. Nevertheless, for a = 1 the radius r is constant and the argument reduces to the angle θ. The JEFs sn, cn and dn are defined as follows: sn(u|m) = y a , (1.66) cn(u|m) = x, (1.67) dn(x|m) = r a , (1.68) where m = k2 is the parameter of the Jacobi elliptic function. In this Thesis the parameter will be denoted as M , to avoid confusion with another parameter that will be labelled m. In the subsequent formulas for the Jacobi elliptic functions the parameter will be implicit. The relationships between the JEFs are: sn2 + cn2 = 1, (1.69) m · sn2 + dn2 = 1. (1.70) When m = 0, the JEFs sn, cn and dn reduce to sin, cos and 1 respectively, whereas when m = 1 the JEFs sn, cn and dn reduce to tanh, sech and sech, respectively. One can then define sc = sncn, ns = 1 sn and so on to obtain a total of 12 Jacobi 24 Table 1.1: Jacobi elliptic functions c0 c2 c4 F M = 0 M = 1 1 1 −(1 +M) M sn sin tanh 2 1−M 2M − 1 −M cn cos sech 3 M − 1 2−M −1 dn 1 sech 4 M −(1 +M) 1 ns cosec coth 5 −M 2M − 1 1−M nc sec cosh 6 −1 2−M M − 1 nd 1 cosh 7 1 2−M 1−M sc tan sinh 8 1−M 2−M 1 cs cot cosech 9 1 −(1 +M) M cd cos 1 10 M −(1 +M) 1 dc sec 1 elliptic functions: sn, cn, dn, ns, nc, nd, sc, cs, sd, ds, cd and dc. The first three JEFs will be called basic JEFs, and the last 9 will be called derived JEFS. The differential equations satisfied by the JEFs are: dsn(u) du = cn(u) · dn(u), (1.71) dcn(u) du = −sn(u) · dn(u), (1.72) dsn(u) du = m · sn(u) · cn(u). (1.73) These properties combined with Eqs. (1.69)-(1.70) give the following form of the differential equation for all 12 Jacobi elliptic functions:( dF dθ )2 = c0 + c2F 2 + c4F 4. (1.74) The coefficients c0, c2 and c4 are given in 1.1 as a function of the elliptic parameter M . The 2 JEFs that are not listed in the table, sd and ds, are not suitable for finding the derived solutions and have hence been omitted. A more extensive list is given in (137) containing 46 different combinations of Jacobi elliptic satisfying an equation of the form given in (1.74), though we will not be using these forms in this Thesis. The forms of the three basic JEFs are given in Fig. (1.1). An incorrect form of Fig. (1.1) was given in (142) and several subsequent papers, that contained 25 formulas for the elliptic modulus, rather than the elliptic parameter, i.e it contained M2 everywhere instead of M . The difference was important only in cases where we combined two different solutions and even in those cases there were only minor quantitative differences with respect to the correct solution, considering that we used values very close to M = 1. All solutions drawn in this Thesis have had this mistake corrected. In Fig. (1.1) we see that the periodicity of each JEF (with the exception of dn at 0) steadily increases asM goes from 0 to 1, diverging to infinity asM approaches 1. The periodicity of JEFs is given as 4K(M), where K is the complete elliptic integral of the first kind. Figure 1.1: Jacobi elliptic functions as a function of u and M : (a) sn, (b) cn, (c) dn, (d) Jacobian amplitude am. 26 1.7 Scope and structure of the Thesis Below I will give a short description of the contents of the various parts of this thesis. In Chapter 2 I present the basics of the F-expansion technique (FET) and the Principle of Harmonic Balance. I then apply the FET to the Nonlinear Schro¨dinger Equation (NLSE) with a third-order nonlinearity to obtain exact spatiotemporal solitary and travelling wave solutions in (3+1) dimensions ((3+1)D). By modifying the method I also obtain similar solutions for the normal dispersion. Finally, I modify the method to obtain solutions for the NLSE with higher-order nonlinearities. In Chapter 3 I apply the FET to the Gross-Pitaevskii (GP) equation. I find that the application of this method depends on solving a Ricatti equation. For sev- eral important systems with a solvable Ricatti equation I obtain numerous solutions. I also describe the methods to obtain novel systems for which the Ricatti equation can be solved. In Chapter 4 I examine the effect of a linear spatial potential on the NLSE. Solutions are obtained for various situations of interest, namely for constant and sinusoidal diffraction and potential as a function of time. In Chapter 5 I explore the application of the FET to new physical systems. Of particular interest are two-component systems such as Manakov systems. I present my results for the case of co-propagating (CO) and counter-propagating (CP) waves. In Chapter 6 I examine the stability of the aforementioned solutions. I apply the standard linear perturbation theory on the NLSE solutions and the GP solu- tions to obtain stability thresholds. The main conclusion is that using dispersion management one can obtain modulationally stable solutions to the NLSE and GP equations. In Chapter 7 I give the main conclusions of my work. 27 1.8 Published papers This section contains a list of all the published research this Thesis was based on, as well as an additional paper still in the process of being finalized: Chapter 2 is based on the following papers: • Wei-Ping Zhong, Rui-Hua Xie, Milivoj Belic´, Nikola Z. Petrovic´, Goong Chen and Lin Yi, ”Exact spatial soliton solutions of the two-dimensional generalized nonlinear Schro¨dinger equation with distributed coefficients,” Phys. Rev. A 78, 023821 (2008), Reference (143). • Milivoj Belic´, Nikola Z. Petrovic´, Wei-Ping Zhong, Rui-Hua Xie and Goong Chen, ”Analytical light bullet solutions to the generalized (3+1)-dimensional nonlinear Schro¨dinger equation,” Phys. Rev. Lett. 101, 0123904 (2008), Reference (142). • Nikola Z. Petrovic´, Milivoj Belic´, Wei-Ping Zhong, Rui-Hua Xie and Goong Chen, ”Exact spatiotemporal wave and soliton solutions to the generalized (3+1)-dimensional nonlinear Schro¨dinger equation for both normal and anoma- lous dispersion,” Opt. Lett. 34, 1609 (2009), Reference (144). • Nikola Z. Petrovic´, Milivoj Belic´ and Wei-Ping Zhong, ”Exact traveling wave and spatiotemporal soliton solutions to the generalized (3+1)-dimensional Schro¨dinger equation with polynomial nonlinearity of arbitrary order,” Phys. Rev. E 83, 026604 (2011), Reference (111). Chapter 3 is based on the following papers: • Nikola Z. Petrovic´, Milivoj Belic´ and Wei-Ping Zhong, ”Exact spatiotempo- ral wave and soliton solutions to the generalized (3+1)-dimensional Gross- Pitaevskii equation,” Phys. Rev. E 81, 016610 (2010), Reference (145). • Nikola Z. Petrovic´, Najdan Aleksic´, Anas Al Bastami and Milivoj Belic´, ”Ana- lytical traveling-wave and solitary solutions to the generalized Gross-Pitaevskii equation with sinusoidal time-varying diffraction and potential”, Phys. Rev. E 83, 036609 (2011), Reference (146). 28 • Anas Al Bastami, Nikola Z. Petrovic´, and Milivoj Belic´, ”Special solutions of the Riccati equation with applications to the Gross-Pitaevskii nonlinear PDE,” Electron. J. Diff. Eqs., Vol. 2010, No. 66, 1 (2010), Reference (147). • Anas Al Bastami, Milivoj Belic´, Danijela Milovic´ and Nikola Z. Petrovic´, ”An- alytical chirped solutions to the (3 + 1)-dimensional Gross-Pitaevskii equation for various diffraction and potential functions,” Phys. Rev. E 84, 016606 (2011), Reference (148). Chapter 4 is based on the following paper: • Nikola Z. Petrovic´, Hussein Zahreddine and Milivoj Belic´, ”Exact spatiotempo- ral wave and soliton solutions to the generalized (3+1)-dimensional nonlinear Schro¨dinger equation with linear potential”, Phys. Scr. 83, 065001 (2011), Reference (149). Chapter 5 is based on the following paper: • Nikola Z. Petrovic´ and Hussein Zahreddine, ”Exact traveling wave solutions to coupled generalized nonlinear Schro¨dinger equations”, Phys. Scr. T149, 014039 (2012), Reference (150). Chapter 6 is based on Ref. (146) and the following paper to be published: • Najdan Aleksic´, Nikola Z. Petrovic´ and Milivoj Belic´, ”General stability anal- ysis of the exact soliton solutions to the generalized nonlinear Schro¨dinger equation and the Gross-Pitaevskii equation”, to be published. 29 Chapter 2 F-expansion technique and its application to the general third-order NLSE 2.1 Introduction The F-expansion technique (FET) is a very well know and common technique for finding solutions to various partial differential equations (PDEs). The main idea is to express the solution of your equation in terms of a polynomial of a certain well- chosen function F . Since we typically deal with nonlinear (NL) partial differential equations of the second order, the well chosen functions will be the solution to a NL differential equation of the form: ( dF dθ )2 = N∑ i=0 ciF i, (2.1) where θ is a parameter that relates to the original variables of the PDE, typically x, y, z and t. One also obtains using chain derivation the following important relationship for F : d2F dθ2 = N∑ i=0 i 2 ciF i−1. (2.2) A famous class of functions often used in solving PDEs are the Weierstrass functions, for which N = 3 and c2 = 0. 30 We will for the most part restrict ourselves in the Thesis to the use of the Jacobi Elliptic Functions (JEFs) for which N = 4 and c1 = c3 = 0. The remaining parameters c0, c2 and c4 depend on the parameter M (in standard literature labeled m) as listed in (1.1). A relation between c0, c2, and c4 holds: c0 + c4 = ±c2. (2.3) Once the form of F is plugged into the solution one obtains a polynomial equa- tion of the form p(F ) = 0, where p is some polynomial function, such that coefficients next to each degree of F can be set to 0. For each degree of F we thus obtain a series of algebraic or ordinary differential equations (ODEs). This is known as the Principle of Homogeneous Balance (PHB). If the form of the solution is well chosen, we obtain novel solutions to the PDE we started from. If the form of the solution was poorly chosen, we obtain mutually contradictory equations or at best the trivial solutions. The system to which we will apply the FET is the general Nonlinear Schro¨dinger equation (NLSE). The general form of the NLSE is: i∂zu+ β(z) 2 N∑ i=1 si∂ 2 xi u+ V (x1, . . . , xn)u+ χ(z)|u|2u = iγ(z)u, (2.4) where u is the wave function, z is the longitudinal coordinate (which can be either spatial or temporal), N is the number of transverse dimensions, xi are the trans- verse coordinates, si = ±1, V is the external potential, which we will usually take to depend only on the transverse variables, χ is the strength of the Kerr nonlinearity and γ is the term for signal gain if positive and signal loss if negative. Here, si = 1 for spatial coordinates and for the temporal coordinate t in the case of the anoma- lous dispersion, while si = −1 for the temporal coordinate in the case of normal dispersion. We will denote in this section (x1, . . . , xN), the set of all transverse vari- ables, as x. This system describes evolution of a slowly-varying wavepacket envelope u(z, x) in a diffractive nonlinear Kerr medium in the paraxial approximation for a wide range of circumstances. The JEFs are suitable functions to utilize in finding exact solutions to the NLSE because from (2.2) one can conclude that the second derivative of the JEF is related to the cube of the original function. When the coefficients are constant, the behavior of solutions to NLSE strongly 31 depends on the dimensionality of the problem. In (1+1)D, as mentioned, one can observe stable localized wavepackets. However, in (2+1)D for the self-focusing non- linearity all localized solutions either spread out with propagation (for input powers less than a critical value) or collapse at a finite distance (for powers above the critical value) (48). This behavior is an example of weak collapse. In (3+1)D one observes the strong collapse: wavepackets collapse at any power, i.e. no power threshold exists. The NLSE is one of the most useful generic mathematical models (151) in many fields of physics. The major scientific interest in this model comes from the fact that in various systems in nonlinear optics and related fields these equations appear naturally. Major interest in the NLSE is piqued by the discovery of solitary wave solutions (6; 4). Stable exact soliton solutions to the NLSE are known only in (1+1)-dimensions ((1+1)D), for the simple reason that the inverse scattering theory (6), responsible for the existence and stability of 1D solitons, only works in (1+1)D. There are no known exact stable solitons in (2+1)D or (3+1)D. Recently, great interest has been generated when it was suggested that the (2+1)D generalized NLSE with varying coefficients may lead to stable 2D solitons (80). The stabilizing mechanism has been the sign-alternating Kerr nonlinearity in a layered medium. A vigorous search for the stabilized periodic solutions of (2+1)D NLSE has been launched (152; 153; 154; 155), however, out of necessity, it has been numerical. In this Thesis, analytical periodic traveling wave and soliton solutions to the NLSE in (3+1)D and several other systems will be presented. 2.2 F-expansion technique We will now apply the general FET (124; 136) and the PHB (119; 123) to find analytical periodic traveling wave solutions to (2.4), for V = 0. We define the complex field u of (2.4) in terms of its amplitude and phase (156): u(z, x) = A(z, x) exp (iB(z, x)) , (2.5) where A and B are real functions. Plugging this form of u into (2.4), cancelling out eiB, and taking the real and imaginary part of the equation separately, we find the 32 following coupled equations for A and B: ∂zA+ 1 2 β ( N∑ i=1 si ( 2∂xiA∂xiB + A∂ 2 xi B )) = γA (2.6) −A∂zB + 1 2 β ( N∑ i=1 si ( ∂2xiA− A (∂xiB)2 )) + χA3 = 0, (2.7) where we have adopted the convention of using indices for partial derivatives, i.e. ∂z ≡ ∂∂z and ∂xi ≡ ∂∂xi . We apply the balance principle with modifications to account for the opposite sign of the time derivative. We seek the traveling wave solutions to (2.6) and (2.7), and assume the functions to be of the form: A(z, x) = f1(z)F (θ) + f2(z)F −1(θ), (2.8) θ = N∑ i=1 ki(z)xi + ω(z), (2.9) B(z, x) = a(z) ( N∑ i=1 six 2 i ) + b(z) ( N∑ i=1 xi ) + e(z), (2.10) where f1, f2, ki, ω, a, b, e are parameter functions to be determined, and F is a JEF satisfying (2.2). The function a is known as the chirp function. The name comes from the applications involving frequency domain. In the spatial domain the chirp functions relates to the curvature of the wave front. The choice of signs si next to a was an important breakthrough that allowed the treatment of normal dispersion within this framework (144). This will be discussed in greater detain in Section (2.3.4). Plugging in Eqs.(2.8)-(2.10) into Eqs.(2.6)-(2.7) and applying d2θF = c2F + 2c4F 3, derived from (2.2), and the corresponding formula for F−1, d2θF −1 = c2F−1+ 2c0F −3, we obtain two complex equations which are polynomial (including negative degrees) with respect to F , dθF = √ c0 + c2F 2 + c4F 4 and the transverse variables xi. One should also note that dθF −1 = −dθF/F 2. Requiring that the equation must hold for arbitrary z, x1, . . . , xN it follows that the equations must also hold for arbitrary F , since the JEFs cover a wide range of parameters. Therefore, the ex- pressions next to each term of the form F j(dθF ) j0xj11 . . . x jN N , where j is an arbitrary integer and ji, i = 0, 1, . . . , N , are arbitrary non-negative integers, all have to be equal to 0. We thus obtain a system of ordinary differential equations (ODEs) and 33 algebraic equations: dfj dz +Naβfj − γfj = 0, (2.11) dki dz + 2kaβ = 0, (2.12) dω dz + β ( N∑ i=1 siki ) b = 0, (2.13) da dz + 2βa2 = 0, (2.14) db dz + 2βab = 0, (2.15) de dz + β 2 (( N∑ i=1 si ) b2 − ( N∑ i=1 sik 2 i ) c2 ) − 3χf1f2 = 0, (2.16) f1 ( β ( N∑ i=1 sik 2 i ) c4 + χf 2 1 ) = 0, (2.17) f2 ( β ( N∑ i=1 sik 2 i ) c0 + χf 2 2 ) = 0. (2.18) We consider the most generic case, in which f1 and f2 are assumed non-zero and β(z) and γ(z) are arbitrary. An alternative case, one where β(z) and χ(z) are arbitrary, but γ(z) is not, will be considered in Sec. 2.3.1.2. Given arbitrary β(z) and γ(z), the following set of exact solutions is found: f1(z) = (α) N/2f10 exp (∫ z 0 γdz ) , (2.19) f2(z) = ǫ √ c0 c4 f1, (2.20) ki(z) = αki0, (2.21) ω(z) = ω0 − α ( N∑ i=1 siki0 ) b0 ∫ z 0 βdz, (2.22) a(z) = αa0, (2.23) b(z) = αb0, (2.24) e(z) = e0 + α 2 (( N∑ i=1 sik 2 i0 ) (c2 − 6ǫ√c0c4)− ( N∑ i=1 si ) b20 )∫ z 0 βdz,(2.25) where ǫ = ±1 and: α(z) = ( 1 + 2a0 ∫ z 0 βdz )−1 (2.26) is the normalized chirp function. The subscript ’0’ denotes the value of the given function at z = 0. 34 One should note the universal influence of the chirp function a on the solutions. In the case when there is no chirp, a0 = 0 and α = 1, the parameters ki and b are all constant. In the presence of chirp they all acquire the prescribed z dependence. The chirp also influences the form of the amplitude A through the dependence of f1, f2 and θ on α. Note that the chirp function figures in the amplitude raised to power N/2, which we term the degree of the chirp. It should also be noted that χ is not arbitrary, but depends on α, β, and γ: χ(z) = −βc4α2−N ( N∑ i=1 sik 2 i0 ) f−20 e −2 ∫ z 0 γdz. (2.27) Hence, in a lossy medium for the nonlinearity coefficient χ will grow exponentially. For N = 2 and no gain or loss (γ = 0), χ stays constant with or without the chirp, while for N 6= 2 the presence of chirp affects χ even without loss or gain. Incorporating (2.19)-(2.25) back into (2.5), we obtain the general periodic trav- eling wave solutions to the generalized NLSE: u = (α)N/2f0e ∫ z 0 γdz ( F (θ) + ǫ √ c0 c4 1 F (θ) ) · (2.28) exp ( i ( a(z) ( N∑ i=1 six 2 i ) + b(z) ( N∑ i=1 xi ) + e(z) )) , where: θ = ω0 + ( N∑ i=1 sikixi ) − ( N∑ i=1 siki ) b0 ∫ z 0 βdz. (2.29) Apart from the solutions given in Eqs. (2.19)-(2.25) one can alternatively assume that f2 = 0, in which case one obtains the exact same equations to which Eqs. (2.19)-(2.25) would reduce for ǫ = 0. Thus, the parameter ǫ in Eq. (2.20) can assume three values: ±1 and 0. As long as one chooses the constants according to the relations listed in Table (1.1) and substitutes the appropriate F (θ) into Eq. (2.8), one obtains exact periodic traveling wave solutions to the generalized (3+1)D NLSE. The parameter M varies between 0 and 1. When M → 0, JEFs degenerate into trigonometric functions, i.e. sn(θ)→ sin(θ), cn(θ)→ cos(θ), dn(θ)→ 1, etc., and the periodic traveling wave so- lutions become the periodic trigonometric solutions. WhenM → 1, JEFs degenerate into hyperbolic functions, i.e. sn(θ) → tanh(θ), cn(θ) → sech(θ), dn(θ) → sech(θ), etc., and the periodic traveling wave solutions become solitary wave solutions. As 35 long as 0 < M < 1 there is no problem with the periodic solutions; one can choose any of the listed functions. However, when M = 0 or M = 1 only some of the functions may be utilized, because of the developing singularities. 2.3 Application of the F-expansion technique We now proceed to apply the results of the previous section to several problems relevant to nonlinear optics. Many of the results for specific equations emerge as special cases of the general equation. We will first study the (3+1)D NLSE with anomalous dispersion to ascertain the general behavior of our solutions. We will then look at the special cases of the (1+1)D and (2+1)D NLSEs. The solutions applied to the (1+1)D NLSE have already been proposed in earlier work so we will only briefly comment on them, while for the (2+1)D NLSEs we will only cover solutions that have no straightforward analogue to the (3+1)D case. We will then consider the case of the (3+1)D NLSE with normal dispersion. 2.3.1 Results for (3+1)-D NLSE with anomalous dispersion The first example of the applications of the FET to the (3+1)D NLSE were published in (142). In this paper the solutions were ambitiously identified as “light bullets”, despite having an obvious extention along the ∑N i=1 siki0xi axis. Alternative methods such as the use of Hirota’s method, or a different form of θ have the potential of producing true light bullet solutions. In this section we cover the generalized NLSE in (3+1)D with distributed coef- ficients (157): i∂zu+ β(z) 2 (∆⊥u+ ∂2t u) + χ(z)|u|2u = iγ(z)u, (2.30) which describes evolution of a slowly-varying wavepacket envelope u(x, y, z, t) in a diffractive nonlinear Kerr medium with anomalous dispersion, in the paraxial ap- proximation. Here z is the propagation (i.e. longitudinal) coordinate, ∆⊥ = ∂2x+∂ 2 y represents the transverse Laplacian, and t is the reduced time. The functions β, χ, and γ stand for the diffraction/dispersion, nonlinearity, and gain coefficients, re- spectively. All coordinates are made dimensionless by the choice of coefficients. The 36 generalized NLSE is of considerable importance, as it describes the full spatiotem- poral optical solitons, called light bullets, in (3+1)D. It is evident that Eq. (2.30) is a special case of Eq. (2.4) for V = 0, N = 3, x1 ≡ x, x2 ≡ y, x3 ≡ t and s1 = s2 = s3 = 1. Thus we take the following form for the solution u that describes traveling waves: u(z, x, y, t) = A(z, x, y, t) exp (iB(z, x, y, t)), (2.31) where: A = f1(z)F (θ) + f2(z)F −1(θ), (2.32) θ = k(z)x+ l(z)y +m(z)t+ ω(z), (2.33) B = a(z)(x2 + y2 + t2) + b(z)(x+ y + t) + e(z). (2.34) Here we have denoted k1 ≡ k, k2 ≡ l and k3 ≡ m. We consider the most generic case, in which f1 and f2 are assumed non-zero. When we apply the FET as described in 2.2 we obtain the following coupled equa- tions: ∂zA+ 1 2 β ( 2 (∂xA∂xB + ∂yA∂yB + ∂tA∂tB) +A(∂2x + ∂ 2 y + ∂ 2 t )B ) = γA, (2.35) −A∂zB + 1 2 β ( (∂2x + ∂ 2 y + ∂ 2 t )A −A ((∂xB)2 + (∂yB)2 + (∂tB)2))+ χA3 = 0. (2.36) Then, plugging in the forms for A and B given in Eqs. (2.32)-(2.34) we obtain: dfj dz + 3aβfj − γfj = 0, (2.37) dki dz + 2kaβ = 0, (2.38) dω dz + β (k + l +m) b = 0, (2.39) da dz + 2βa2 = 0, (2.40) db dz + 2βab = 0, (2.41) de dz + β 2 ( 3b2 − (k2 + l2 +m2) c2)− 3χf1f2 = 0, (2.42) f1 ( β ( k2 + l2 +m2 ) c4 + χf 2 1 ) = 0, (2.43) f2 ( β ( k2 + l2 +m2 ) c0 + χf 2 2 ) = 0. (2.44) 37 We will now proceed to solve Eq. (2.30) for two cases. In the first case we assume that β(z) and γ(z) are arbitrary functions and χ(z) is ultimately determined to have a specific form as an integrability condition. In the second case, β(z) and χ(z) are arbitrary and γ(z) is determined. The case of γ(z) and χ(z) being arbitrary has not been solved, due to the dependence of most equations on β(z). 2.3.1.1 Case 1: Arbitrary β(z) and γ(z) We will first consider the case where β(z) and γ(z) are arbitrary. In this case, we obtain the solutions to Eqs. (2.37)-(2.44) as follows: f1(z) = (α) 3/2f0 exp (∫ z 0 γdz ) , (2.45) f2(z) = ǫ √ c0 c4 f1, (2.46) k(z) = αk0, (2.47) l(z) = αl0, (2.48) m(z) = αm0, (2.49) ω(z) = ω0 − α(k0 + l0 +m0)b0 ∫ z 0 βdz, (2.50) a(z) = αa0, (2.51) b(z) = αb0, (2.52) e(z) = e0 + α 2 ( (k20 + l 2 0 +m 2 0)(c2 − 6ǫ √ c0c4)− 3b20 ) ∫ z 0 βdz, (2.53) where ǫ = ±1 and: α(z) = ( 1 + 2a0 ∫ z 0 βdz )−1 (2.54) is the normalized chirp function. The subscript ’0’ denotes the value of the given function at z = 0 and f0 is defined to be f10. The Jacobi functions one may use in this solution, as well as the corresponding c0, c2 and c4 are given in Table (1.1). In order to avoid singularities when β is the cosine function, we must have the condition on the initial value of the chirp |a0| < 1/2. The integrability condition for χ reduces to: χ = −βc4(k20 + l20 +m20)f−20 e−2 ∫ z 0 γdz/α. (2.55) As in Sec. 2.2, one obtains assuming f2 = 0 the exact same equations to which Eqs. (2.45)-(2.53) would reduce for ǫ = 0. Thus, the parameter ǫ in Eq. (2.20) can assume three values: ±1 and 0. 38 Incorporating Eqs. (2.45)-(2.53) solutions back into (2.31), we obtain the general periodic traveling wave solutions to the generalized NLSE: u = (α)3/2f0e ∫ z 0 γdz ( F (θ) + ǫ √ c0 c4 1 F (θ) ) · exp i (a(x2 + y2 + t2) + b(x+ y + t) + e) , (2.56) where: θ = ω0 + kx+ ly +mt− (k + l +m)b0 ∫ z 0 βdz. (2.57) We should note that for M = 1 the solution described by Eq. (2.56) describes spatially extended spatio-temporal (ST) solitons. Even though the amplitude A as a function of the transverse variable θ is localized, it is not when viewed in the plane of transverse coordinates x and y. This is easily seen if one rotates the x and y axes about the z axis for some angle α, to arrive at a set of new coordinates x′ and y′. By choosing the angle as tan(α) = −k/l, the variable θ will not contain y′, and by choosing tan(α) = k/l, it will not contain x′. Thus the amplitude A will not explicitly depend on y′ (or x′) and the soliton will be extended along the y′ (or x′) axis. Hence, the solutions obtained with the present method cannot be of the light bullet type (142). As an example, we present some of the periodic wave and light bullet soliton solu- tions, taking the diffraction/dispersion coefficient β to be of the form β = β0 cos (kbz) and the gain/loss coefficient γ to be a small constant. This choice leads to alternat- ing regions of positive and negative values of both β and χ, which is required for an eventual stability of soliton solutions. We will typically take kb = 1. In Figs. (2.1)-(2.2), we depict the periodic wave solutions made up from the single F functions, F = sn and F = cn, functions 1 and 2 from Table (1.1). In Figure (2.1) there is no chirp (a0 = 0), while in Figure (2.2) there is chirp (a0 = 0.1). For both figures ǫ = 0. We see that in the absence of chirp the function is periodic in the transverse direction, while in the case of the chirp there is stretching in the transverse direction. Also, the amplitude of the function is affected by the chirp. These effects grow stronger with the increase of chirp. An important feature that distinguishes our solutions from the others reported in the literature (157; 140), apart from the dimensionality, is this appearance of the general spatiotemporal chirp function in both phase and amplitude. 39 Figure 2.1: Periodic traveling wave solutions as functions of the propagation dis- tance, for a0 = 0 (without chirp) and ǫ = 0. Intensity |u|2 for: (a) F = sn and (b) F = cn, presented as functions of k0x+ l0y+m0t and z. Coefficients: β(z) = cos(z), γ(z) = γ0 = 0, M = 0.9999, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0. Figure 2.2: Periodic traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.1), except for a0 = 0.1. 40 Figure 2.3: Combined intensity distributions of the periodic waves 1 and 4 from Table (1.1) (F = sn and F = ns), as functions of the propagation distance, with ǫ = 1 for: (a) a0 = 0 (no chirp) and (b) a0 = 0.1 (with chirp). Other parameters are the same as in Fig. (2.1). Figure (2.3) shows the periodic wave solution made from the combination of the F functions 1 and 4 for ǫ = 1, also without and with the chirp. As can be seen, the presence of ǫ significantly changes the nature of solutions. The parameterM significantly affects the periodicity of the solutions. In Figure (2.4) we see the effect of parameter M on the solutions. We can clearly see that as M decreases, so does the periodicity of our travelling waves. It should be noted that the Jacobi elliptic function also changes shape with the change of M . On the other hand, as M approaches 1 the period of the JEF stretches to infinity and one obtains solitary waves. Figures (2.5)-(2.7) repeat the same sequence of plots as Figs. (2.1)-(2.3), but show the light bullet soliton solutions instead. The soliton solution is similar to a single period of the periodic wave solutions for M close to 1. The obtained solutions are dependent on a number of parameters and it is instructive to see how the form of the solution changes with the change of various parameters. The key parameter is the form of the diffraction/dispersion coefficient β. When β is constant the chirp function α converges to 0 or to infinity. We see in Fig. (2.8) that based on the signs of a0 and β, that the chirp will either attenuate or dissipate the solution. Also, these solutions tend to be unstable. This is why recent research has focused on periodic forms of β, and as a consequence of (2.55) also of 41 Figure 2.4: Periodic traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.1)(a), except for: (a) M = 0.999, (b) M = 0.99, (c) M = 0.9 and (d) M = 0.5. Figure 2.5: Solitary wave solutions without chirp. The setup and parameters are as in Fig. (2.1), except for M = 1. 42 Figure 2.6: Solitary wave solutions with chirp. Setup is the same as in Fig. (2.5), except for a0 = 0.1. Figure 2.7: Combined intensity distributions of the solitary wave solitons 1 and 4 (F = sn and F = ns). The setup is as in Fig. (2.3) except for M = 1. 43 Figure 2.8: Solitary traveling wave solutions, as functions of the propagation dis- tance. The setup and parameters are the same as in Fig. (2.5), except for F = cn and: (a) β(t) = 1, b0 = 0.1, a0 = 0, (b) β(t) = 1, b0 = 0.1, a0 = 0.02 (c) β(t) = −1, b0 = 0.1, a0 = 0 and (d) β(t) = −1, b0 = 0.1, a0 = 0.02. χ. Making β and χ periodic is known as the dispersion management. In Fig. (2.9), the graph of the integrability condition is given for a sinusoidal β. We see that for N = 1 and N = 3 the form of χ is deformed due to the presence of chirp. However, the deformation is not significant enough to produce practical difficulties in constructing χ. For N = 2, the integrability condition is not affected by the chirp and the graph of χ remains sinusoidal for all values of a0. In Fig. (2.10) the effect of the chirp is shown. It is evident that as the chirp increases the solitary wave increases on one side and decreases on the other. The two sides are determined by the sign of the initial value of chirp a0. In Fig. (2.11) we see the effect of parameter b0 on the form of the solutions. 44 Figure 2.9: The nonlinearity streingth given the integrability condition, as a function of the the propagation distance z and the initial chirp a0. The parameters are the same as in Fig. (2.5), except for F = cn and: (a) N = 1 and (b) N = 3. Roughly, the parameter b determines the intensity with which β will be expressed in the final form of the solutions. As can be seen from Eq. (2.50), when b0 = 0 ω will be constant regardless of the form of β and will ultimately lead to the solitary wave moving in a straight line. We see from Fig. (2.11)(b) that the effect of chirp, i.e. the change of amplitude is maintained. In Figure (2.12), we see the effects of the presence of loss in the system γ = −0.05. As we can see, the solutions roughly maintain their shape while their am- plitude decays. We also see that solutions other than 1 and 2 entail singularities, hence they are not of as much interest to us. In Figure (2.13), we see the same effects for travelling wave solutions. The singularities occur in two sets of alternating regions, corresponding to whether sn or cn is in the denominator. In Figure (2.14) we see the effects of chirp on decaying functions for F = sn and F = cn in Fig. (2.12) and Fig. (2.13). Finally, in Fig. (2.15) we see the graphs of the phase of our solutions, with and without chirp. 45 Figure 2.10: Solitary traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.5), except for: (a) F = sn, a0 = 0.2 (b) F = cn, a0 = 0.2, (c) F = sn, a0 = −0.1 and (d) F = cn, a0 = −0.1. 46 Figure 2.11: Solitary traveling wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.5), except for F = cn and: (a) b0 = 0 (b) b0 = 0, a0 = 0.1, (c) b0 = 0.5 and (d) b0 = −1. 47 Figure 2.12: Solitary wave solutions without chirp. The intensities |u|2 are plotted of solutions: (a) 1 (F = sn), (b) 2 (F = cn), (c) 7 (F = sc), (d) 8 (F = cs), (e) 4 (F = ns) and (f) 5 (F = nc). The parameters are otherwise the same as in Fig. (2.1), except M = 1 and γ(z) = −0.05. 48 Figure 2.13: Traveling wave soliton solutions without chirp. The parameters are the same as in Fig. (2.12), except M = 0.9999. 49 Figure 2.14: Solitary and travelling wave soliton solutions with chirp. The parame- ters in (a) and (b) are the same as in Fig. 2.12 (a) and (b), and the parameters in (c) and (d) are the same as in 2.13 (a) and (b), except a0 = 0.1. 50 Figure 2.15: Phase of the solutions B as functions of the propagation distance and one transverse variable, assuming y = t = 0. The parameters are the same as in Fig. (2.5), except for F = cn and: (a) a0 = 0, and (b) a0 = 0.1. 2.3.1.2 Case 2: Arbitrary β(z) and χ(z) We now assume that β(z) and χ(z) are arbitrary functions and proceed to solve Eqs. (2.37)-(2.44). Parameters f2, a, b, k, l, m, ω and e have the same form, as given in Eqs. (2.46)-(2.53). Parameters f1, and γ have the following values: f1 = α √ β(z) χ(z) (k20 + l 2 0 +m 2 0)c4, (2.58) γ(z) = 3a(z)β(z) + df1 dz f1 , (2.59) where α is the chirp function given in Eq. (2.54). Given that c4 = M for F = sn and c4 = −M for F = cn, it follows that β(z)/χ(z) must be positive in the case of dark solitons and negative in the case of bright solitons. This set of solutions is significant in that, given adequate gain or loss, one can produce a wide range of solutions. For constant or proportional β and χ one obtains solutions similar to those in Sec. 2.3.1.1 except that the degree of chirp is 1 instead of 3/2. Still, the two sets of solutions nevertheless appear very similar qualitatively. The integer degree of chirp makes even solutions with negative chirp possible, which will be discussed in Sec. 2.3.3. It is, however, possible to introduce sinusoidal modifications in either β or χ giving us an oscillation of amplitude even in the complete absence of chirp. In Fig. 51 Figure 2.16: Solitary wave solutions for β(z) = 1 + sin(z), f0 = 1, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. The remaining parameters are: (a) a0 = 0, F = sn and χ(z) = 1, (b) a0 = 0, F = cn and χ(z) = −1, (c) a0 = 0.02, F = sn and χ(z) = 1, and (d) a0 = 0.02, F = cn and χ(z) = −1. (2.16) we see the effect of a shifted sinusoidal value of β for a constant χ. The solitary wave keeps travelling to one side, due to the function γ required to maintain this motion being similar to the cot function, and it doesn’t vary too much as a function of chirp, as shown in plot (a) of Fig. (2.18). The travelling waves are not shown, but they exhibit the same trends as solutions in Sec. 2.3.1.1: waves travel in parallel without chirp and spread out with chirp. In Fig. (2.17) we see the effects of varying χ while keeping β constant. In order to avoid singularities, we took χ(z) = 2 + sin(z) for the dark soliton and χ(z) = −(2 + sin(z)) for the bright soliton in our research. The key difference 52 Figure 2.17: Solitary wave solutions for β(z) = 1, χ(z) = −(2+ sin(z)) and F = cn. Other parameters are the same as in Fig. (2.16), except for: (a) a0 = 0 and (b) a0 = 0.02. between the solitons in Fig. (2.17) and those in Fig. (2.16) is that solitons in Fig. (2.17) travel in perfectly straight lines. We show only the bright solitary waves, with and without chirp. Finally, in Fig. (2.18) we show the values of γ needed to achieve both regimes. For the first case (Fig. (2.16)) the value of γ is sharply discontinuous at points where f2 reached 0, whereas for the second case (Fig. (2.17)) there is a periodic dependance. The values of γ are given for the cases of a0 = 0 and a0 = 1 and not much difference is seen at higher values of z in both cases. 2.3.2 Results for (1+1)D NLSE We will now cover a few other special cases of the general system given in Sec. 2.2. The formulas for the (1+1)D equation can be obtained from Equations (2.19)-(2.25) by taking N = 1 and s1 = 1. Note that the degree of the chirp function is 1/2 in the expression for f1 and f2 given in Eq. (2.19). It turns out that Kruglov et al. first found solutions to the (1+1)D NLSE that use the Jacobi elliptic functions in Ref. (140). In this paper, the authors postulated the existence of the quadratic term in the phase of the solution, though they assumed that the linear term, i.e. b(z) in our notation, was 0. The authors recognized the 53 Figure 2.18: Values of γ(z) needed to achieve solution, as given in Eq. (2.59). Parameters in (a) are the same as given in Fig. (2.16), whereas the parameters in (b) are the same as given in Fig. (2.17). In both graphs a0 = 0 for the lower plot and a0 = 1 for the upper plot. significance of the Jacobi elliptic functions in finding solutions to the (1+1)D NLSE and they list several solutions that utilize the different JEFs. The solutions given in (140) match our solutions for the (1+1)D NLSE, including the degree of chirp. In the appendix, these authors use the so-called autonomous principle to argue why it is not possible to include phase coefficients of degree higher than 2 with respect to the transverse variables for the (1+1)D NLSE with Kerr nonlinearity. 2.3.3 Results for (2+1)D NLSE The (2+1)D case was extensively covered in (143), which preceded the publication of (142) and the other papers covered in this Thesis. We will summarize here all the results not covered in Sec. 2.3.1. The formulas for the (2+1)D equation can be obtained from Eqs. (2.19)-(2.25) by taking N = 2 and s1 = s2 = 1. In (143), the formula for the chirp as a function of β(z) was not explicitly stated. Instead, the chirped and the unchirped cases were treated separately, wherein for the chirped case the function a(z) was simply postulated, and the formula for β was derived via the relationship: β(z) = − 1 2a2 da dz . (2.60) The order of the chirp function in Eq. (2.19) is 1. Since the formula doesn’t contain 54 the square root, it is possible for the chirp function to be negative, unlike in the cases of the (1+1)D and (3+1)D NLSE. If the chirp function crosses the x-axis it is evident from Eq. (2.60) that there is a singularity in β. However, the only parameters in which β appears, ω and e, can be rewritten as follows: ω(z) = ω0 − (k0 + l0)b0q, (2.61) e(z) = e0 + 1 2 ( (k20 + l 2 0 +m 2 0)(c2 − 6ǫ √ c0c4)− 3b20 ) q, (2.62) where from Eq. (2.26) one obtains: q ≡ α ∫ z 0 β(z)dz = 1− α(z) 2a0 = 1− a(z) a0 2a0 . (2.63) This parameter q that is common to ω and e will prove important in Chapter 3. In Fig. (2.19) we see novel solution obtained when the chirp is of the sinusoidal form. As expected, the pattern for a solitary-wave solution is repeated and stretched out due to the effect of chirp. Figures (b) and (d) were presented for the first time in (143), although here they are shown on a wider scale. Note that we cannot begin integration at 0, but must do so at a point without singularity, for example z = π/2. Thus the index ’0’ represents the values of the given parameter at z = π/2, and hence a0 will stand in this case for the amplitude of the chirp at that point. In Fig. (2.20) we see the effects of chirp on the amplitude for β(z) = − 1 2(1+z)2 . If a0 = 0 then there is no chirp, but if a0 = 1 the chirp becomes a(z) = 1 + z. Thus, since the square of the chirp factors in the form of the amplitude, in figure (b) we see the quadratic increase of amplitude in both +z and −z directions compared against a sinusoidal form of the gain function γ(z). The figures in this graph were also presented for the first time in (143) albeit over a much smaller range. 2.3.4 Results for (3+1)-D NLSE with normal dispersion We now consider the same NLSE in (3+1)D but with normal dispersion, instead of the anomalous dispersion: i∂zu+ β(z) 2 (∆⊥u− ∂2t u) + χ(z)|u|2u = iγ(z)u. (2.64) All of the quantities in this equation are the same as described in 2.3.1. This equation is a special case of (2.4) for V = 0, N = 3, x1 ≡ x, x2 ≡ y, x3 ≡ t, s1 = s2 = 1 and 55 Figure 2.19: Solutions for a sinusoidal chirp of the (2+1)-D NSLE. The parameters are: α(z) = sin(z), γ(z) = sin(z), f0 = 1, b0 = 1, e0 = 0, k0 = l0 = 1, ω0 = 0, ǫ = 0 and: (a) M = 1, F = sn, (b) M = 1, F = cn, (c) M = 0.999, F = sn and (d) M = 0.999, F = cn. 56 Figure 2.20: Solutions for β(z) = − 1 2(1+z)2 , γ(z) = cos(z)/2, F = sn, M = 1. Other parameters are the same as in Fig. (2.19) except for: (a) α = 0 and (b) α = 1 + z. s3 = −1. The generalization described in Sec. 2.2 that allowed the treatment of the (3+1)-D NLSE was first described in (144). Applying the solution form for u, as in (2.31), we assume the form of A and B as follows: A = f1(z)F (θ) + f2(z)F −1(θ), (2.65) θ = k(z)x+ l(z)y +m(z)t+ ω(z), (2.66) B = a(z)(x2 + y2 − t2) + b(z)(x+ y + t) + e(z), (2.67) where f , g, k, l, m, ω, a, b, e are parameter functions to be determined, and F is a Jacobi elliptic function (JEF). As in (2.3.1) k1 ≡ k, k2 ≡ l and k3 ≡ m. These equations are a special case of Eqs. (2.8)-(2.10) and differ from Eqs. (2.32)-(2.34) only in the sign next to t2 in (2.67). As in (2.3.1), we consider the most generic case, in which f and g are assumed non-zero and β(z) and γ(z) are arbitrary. We also assume k2+l2−m2 6= 0, otherwise the only solution for non-zero χ is f = g = 0. When we apply the FET as described 57 in 2.2 we obtain the following solutions: f1(z) = (α) 3/2f0 exp (∫ z 0 γdz ) , (2.68) f2(z) = ǫ √ c0 c4 f1, (2.69) k(z) = αk0, (2.70) l(z) = αl0, (2.71) m(z) = αm0, (2.72) ω(z) = ω0 − α(k0 + l0 −m0)b0 ∫ z 0 βdz, (2.73) a(z) = αa0, (2.74) b(z) = αb0, (2.75) e(z) = e0 + α 2 ( (k20 + l 2 0 −m20)(c2 − 6ǫ √ c0c4)− b20 ) ∫ z 0 βdz, (2.76) where α = (1 + 2a0 ∫ z 0 βdz)−1 is the normalized chirp function. The subscript ‘0’ denotes the value of the given function at z = 0 and f0 is defined to be f10. The constants c0, c2, and c4 are related to the elliptic modulus of JEFs. A parameter ǫ = ±1 is introduced in Eqs. (2.68)-(2.76), to distinguish the two present possibilities. One should note the universal influence of the chirp function α on the solutions, similar to the influence of chirp in the case of normal dispersion. In the case when there is no chirp, a0 = 0 and α = 1, the parameters k, l, m and b are all constant. In the presence of chirp they all acquire the prescribed z dependence that also influences the form of the amplitude A through the dependence of f , g, and θ on α. As in the case for the anomalous dispersion, χ is not arbitrary, but depends on α, β, and γ: χ = −βc4(k20 + l20 −m20)f−20 e−2 ∫ z 0 γdz/α. (2.77) Incorporating these solutions back into Eqs. (2.65)-(2.67), we obtain the general periodic traveling wave solutions to the generalized NLSE: u = (α)3/2f0e ∫ z 0 γdz ( F (θ) + ǫ √ c0 c4 1 F (θ) ) exp i(a(x2+y2−t2)+b(x+y+t)+e), (2.78) where: θ = ω0 + kx+ ly +mt− (k + l −m)b0 ∫ z 0 βdz. (2.79) As in Sections 2.2 and 2.3.1, assuming f2 = 0, one obtains ǫ = 0. Thus, the 58 parameter ǫ in Eq. (2.78) can assume three values: ±1 and 0. As in Sec. 2.3.1 we can choose the form of the solutions from Table (1.1). We present some of the travelling wave and solitary solutions, taking the diffrac- tion/dispersion coefficient β to be of the form β = β0 cos (kbz) and the gain/loss coefficient γ to be a small negative constant. In Figs. (2.21)-(2.22) we depict the solitary wave solutions made up from the single F functions 1 and 2 from the table, without and with the chirp, for ǫ = 0. For the same values of k0, l0 and m0 the overall profile is narrower in the case of normal dispersion, due to a change of sign in m0, but the overall nature of the solutions is not changed. It is worth noting that for m0 = k0 + l0, the profile of the solitary wave is straight even for nonzero b0, i.e. it resembles Fig. (2.11)(a) and (b). For m0 > k0 + l0 the profile oscillates in the other direction, resembling Fig. (2.11)(d) for the case without chirp. Figure 2.21: Solitary wave solutions as functions of the propagation distance, for a0 = 0 (without chirp) and ǫ = 0. Intensity |u|2 of: (a) solution 1 and (b) solution 2 from Table (1.1), presented as functions of k0x+ l0y+m0t and z. Other parameters are: β(z) = cos(z), γ(z) = γ0 = −0.05, M = 1, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0. Figures (2.23)-(2.24) repeat the same sequence of plots as Figures (2.21)-(2.22), but show the traveling wave solutions instead. The soliton solution is similar to a single period of the periodic wave solutions for M close to 1. The most striking feature of the solutions for the normal dispersion is their 59 Figure 2.22: Solitary wave solutions with the chirp, as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.21), except for a0 = 0.1. Figure 2.23: Traveling wave solutions without chirp. The setup and parameters are as in Fig. (2.21), except for M = 0.9999. 60 Figure 2.24: Traveling wave solutions with chirp. Setup is the same as in Fig. (2.23), except for a0 = 0.1. similarity to those for the anomalous dispersion. With the modifications of only a few parameters, the solutions for this fundamentally different system from that of the NLSE with anomalous dispersion were obtained. 2.4 Adaptation of the F-expansion technique for higher order nonlinearities We now turn to the question of generalizing these results to higher order nonlin- earities, i.e. we present analytical periodic traveling wave and soliton solutions for the Kerr nonlinearity of arbitrary high order. Solitons involving higher orders of nonlinearity were first developed in (158) and finding solutions to the higher order NLSEs, especially the Cubic-Quintic equation and the Septic equation, remains an ongoing area of research. The results in this section were first developed in (111). We are interested in the generalized NLSE in (3+1)D with distributed coeffi- cients (157): i∂zu+ β(z) 2 ∆u+ χ1(z)|u|2u+ · · ·+ χn(z)|u|2nu = iγ(z)u, (2.80) which describes evolution of a slowly-varying envelope u(x, y, t, z) in a diffractive and dispersive nonlinear medium, with an arbitrarily large order of nonlinearity. Here z is the propagation coordinate and ∆ = ∂2x + ∂ 2 y + s∂ 2 t represents the generalized 3D 61 transverse Laplacean, in which x and y are the transverse spatial coordinates, and t is the reduced time. As in the previous systems, all coordinates are made dimensionless by the choice of coefficients and functions β and γ stand for the diffraction/dispersion and gain coefficients, respectively. In this section we assume s = 1, i.e. the case of anomalous dispersion, although the results are easily generalizable to the case of normal dispersion as well. The functions χm for m = 1, 2, . . . , n stand for the nonlinearities of orders up to 2n+1. For n = 1 one has the simple Kerr nonlinearity, for n = 2 the cubic-quintic, for n = 3 the septic, and so on. Apart from the study of cubic-quintic and septic systems, another motivation to look into exact solutions of the generalized NLSE with high-order Kerr nonlinearity comes from the fact that such a nonlinearity is an excellent approximation to the saturable nonlinearity: 1 1 + sI ≈ 1− sI + (sI)2 − (sI)3 + . . . , (2.81) which is an important generic model, but for which there are no known exact solu- tions. Following the standard procedure as described in Sec. 2.2, we write the complex field u of (2.80) in terms of its amplitude and phase: u(z, x, y, t) = A(z, x, y, t) exp (iB(z, x, y, t)) . (2.82) Substituting u into Eq. (2.80), the following coupled equations are obtained: ∂zA+ 1 2 β(2∂xA∂xB + 2∂yA∂yB + 2∂tA∂tB+ A∆B) = γA, (2.83) −A∂zB + 1 2 β ( ∆A− A (∂xB)2 − A (∂yB)2 − A (∂tB)2 ) + χ1A 3 + · · ·+ χnA2n+1 = 0. (2.84) As in previous sections, we apply the balance principle (119; 123) and the F- expansion technique (124; 136), as developed in (156), with modifications to ac- count for the higher order nonlinearities. We seek the traveling wave solutions to Eqs. (2.83)-(2.84), and assume the A and B functions to be of the form: A = f1(z)F 1 n (θ) + f2(z)F − 1 n (θ), (2.85) θ = k(z)x+ l(z)y − Ω(z)t+ φ(z), (2.86) B = a(z)(x2 + y2 + t2) + b(z)(x+ y + t) + e(z), (2.87) 62 where f1, f2, k, l, Ω, φ, a, b, e are parameter functions to be determined, and F is a Jacobi elliptic function (JEF). In this section we will use Ω = −m instead of m in order to obtain travelling waves in the positive direction. These solutions resemble the solutions used in the original F-expansion technique (Eqs. (2.32)-(2.34)), except for the power of the function F , which is now 1/n instead of 1. The change in the power is crucial, allowing the higher order nonlinearities to be accounted for. The power has to be such that the highest-order term from Laplacean matches the highest-order nonlinearity. We substitute Eqs. (2.85)-(2.87) into Eqs. (2.83)-(2.84) and in accordance with the F-expansion technique require that xjF 2p−1 n , yjF 2p−1 n , tjF 2p−1 n , (j = 0, 1, 2, p = −n, . . . , n+1) and √c0 + c2F 2 + c4F 4 of each term be separately equal to zero. The constants c0, c2, and c4 are related to the elliptic modulus M of JEFs as given in Table (1.1). After multiplying the expressions and factoring out common factors of f1 and f2, a system of first-order ordinary differential equations is obtained for the parameter functions: dfj dz + 3aβfj − γfj = 0, (2.88) dk dz + 2kaβ = 0, (2.89) dl dz + 2laβ = 0, (2.90) dΩ dz + 2Ωaβ = 0, (2.91) dφ dz + β(k + l − Ω)b = 0, (2.92) da dz + 2βa2 = 0, (2.93) db dz + 2βab = 0, (2.94) −de dz − 3 2 βb2 + qc2 + n∑ i=1 χi ( 2i+ 1 i+ 1 ) = 0, (2.95) where j = 1, 2, q = β(k 2+l2+Ω2) 2n2 , and by definition χm ≡ χmfm1 fm2 , m = 1, . . . , n. A number of algebraic relations involving χm is also obtained: n∑ i=1 χi ( 2i+ 1 i+ p ) = 0, (2.96) χn−1 + (2n+ 1)χn − (n− 1)qw = 0, (2.97) χn + (n+ 1)qw = 0, (2.98) 63 where w = c0 ( f1 f2 )n = c4 ( f2 f1 )n , and p = 2, . . . , n − 1. Indeed, Equations (2.95), (2.97) and (2.98) are obtained from the terms next to F 1 n , F 2− 1 n and F 2+ 1 n , respec- tively. Note that Eq. (2.97) appears only for n > 1 and Eq. (2.96) only for n > 2. The binomial coefficient ( 2i+1 i+p ) is defined to be 0 for i+ p > 2i+ 1. Equations (2.88)-(2.98) resemble Eqs. (2.11)-(2.18), the corresponding equations in Sec. 2.2 for the special case ofN = 3 and V = 0 , except that the equations for e(z) and χn(z) have now been generalized. By solving these equations self-consistently, one obtains a set of conditions on the coefficients and parameters, necessary for Eq. (2.80) to have exact traveling wave solutions. We consider the most generic case, in which β(z) and γ(z) are arbitrary. We first solve Eqs. (2.88)-(2.94), to obtain expressions for f1, f2, k, l, Ω, φ, a and b. From the condition on w it follows f2 = ǫf1 2n √ c0 c4 and so w = ǫn √ c0c4, where ǫ = ±1. We then proceed to solve for χm recurrently starting from m = n and ending at m = 1. We easily obtain: χm = rmqw (2.99) and: n∑ i=1 χi ( 2i+ 1 i+ 1 ) = rqw, (2.100) where the r-parameters r, r1, . . . , rn are integer functions of n. Although it is difficult to find generic formulas for r, r1, . . . , rn, in principle it is easy to calculate these parameters recursively for any given n: rn = −(n+ 1), (2.101) rn−1 = (n− 1)− (2n+ 1)rn = 2n(n+ 2), (2.102) rm = − n∑ i′=m+1 ri′ ( 2i′ + 1 m+ i′ + 1 ) , m = 1, ..., n− 2, (2.103) r = n∑ i=1 ri ( 2i+ 1 i+ 1 ) . (2.104) The values of coefficients r and r1, . . . , rn are given in Table (2.1). In the end, the 64 Table 2.1: Values of r-parameters n r r1 r2 r3 r4 r5 r6 r7 1 -6 -2 2 18 16 -3 3 -38 -66 30 -4 4 66 192 -156 48 -5 5 -102 -450 570 -300 70 -6 6 146 912 -1659 1312 -510 96 -7 7 -198 -1666 4116 -4536 2590 -798 126 -8 8 258 2816 -9072 13248 -10340 4608 -1176 160 9 -326 -4482 18252 -34056 34650 -20826 7602 -1656 10 402 6800 -34155 79200 -101530 78624 -38325 11840 n r8 r9 r10 8 -9 9 198 -10 10 -2250 240 -11 65 following set of exact solutions is found: f1(z) = (α) 3/2f0e ∫ z 0 γdz, (2.105) f2(z) = ǫ 2n √ c0 c4 f1, (2.106) k(z) = αk0, (2.107) l(z) = αl0, (2.108) Ω(z) = αΩ0, (2.109) φ(z) = φ0 − α(k0 + l0 − Ω0)b0 ∫ z 0 βdz, (2.110) a(z) = αa0, (2.111) b(z) = αb0, (2.112) e(z) = e0 + α ( (k20 + l 2 0 + Ω 2 0) 2n2 (c2 + rǫ n√c0c4)− 3b 2 0 2 )∫ z 0 βdz. (2.113) As for previous systems, α(z) = (1 + 2a0 ∫ z 0 βdz)−1 denotes the normalized chirp function, subscript 0 denotes the value of the given function at z = 0 and f0 = f10. Since A has to be real, for even n we must have F > 0 at all times, so that F 1/n is real. This restricts the range of allowed solutions in terms of JEFs. Another restriction involves the nonlinearity coefficients, which by the solution procedure are found related to β and γ: χm = ǫn+mrmβα 2−3m 2n2f 2m0 2n √ cn+m4 c n−m 0 exp ( −2m ∫ z 0 γdz ) , (2.114) where m = 1, 2, . . . , n. This relation should be understood as an integrability con- dition on Eq. (2.80), similar to the formula for χ obtained in Eq. (2.27) for the system studied in Sec. 2.2. Note that the nonlinearity coefficients χm are directly proportional to the rm parameters, while the only parameter to explicitly appear in the solutions is r, which appears in Eq. (2.113). We consider separately the case f2 = 0. In this case the solutions are given in terms of single JEFs. We then have one negative exponent of F in the term ∆A of Eq. (2.83), namely (n − 1)f1c0F−2+ 1n . However, all other terms in (2.83) will be with a positive degree of F , hence for n 6= 1 we must have c0 = 0. Then one finds χn = −q(n + 1)c4f−2n1 and all other χm = 0, m = 1, . . . , n − 1. The correct expressions are still obtained from Eqs. (2.105)-(2.113), provided one takes ǫ = 0. The case n = 1 was covered in Sec. 2.3.1; the expression for χ Eq. (2.55) remains unchanged, c0 need not be 0, and the correct expressions are again contained in Eqs. 66 (2.105)-(2.113), as long as one takes ǫ = 0. All of these restrictions do not bode well for the application of the present solution method to the saturable Kerr-like nonlinearity, which was one of our original aims. For the time being, this important model remains nonintegrable (159). 2.4.1 Results for n = 1 and n = 2 We now apply the results obtained to a few specific cases. For n = 1 we have r1 = −2 and r = −6. Hence, Eqs. (2.105)-(2.113) reduce to the corresponding equations (2.45)-(2.53) obtained in Sec. 2.3.1. For n = 2, corresponding to the cubic-quintic model, we obtain r1 = 16, r2 = −3 and r = 18. The only soliton solution found so far which satisfies both F > 0 and c0 = 0 for f2 = 0 is the bright soliton solution with M = 1 and F = sech. This solution is presented in Fig. (2.25). Figure 2.25: Soliton solutions for the cubic-quintic model n = 2 as a function of the propagation distance, for: (a) a0 = 0 (without chirp) and (b) a0 = 0.1 (with chirp) for the F = sech solution. Intensity |u|2 presented as a function of k0x + l0y − Ω0t and z. Other parameters are: β(z) = cos(z), γ(z) = γ0 = −0.05, M = 1, b0 = 1, e0 = 0, ǫ = 0, k0 = l0 = −Ω0 = 1, φ0 = 0. 67 2.4.2 Results for n = 3 For n = 3, corresponding to the septic model, we have r1 = −66, r2 = 30, r3 = −4 and r = −38. We obtain the bright soliton solution for f2 = 0, which looks very much like the solution seen in Fig. (2.25). The only noticeable difference is the transverse stretching of the wave. This is due to the fact that the function F 1 n falls less rapidly as the argument decreases for larger n. For the septic model we also obtain solutions for ǫ = 1. These correspond to the dark solitons, with F = tanh. Since c0 = 0 is no longer required, one can find both solitary (M = 1) and periodic (M < 1) traveling wave solutions. These solutions are shown in Figs. (2.26)-(2.27). Figure 2.26: Soliton solutions for the septic model as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.25), except for n = 3, F = tanh and ǫ = 1. 68 Figure 2.27: Periodic traveling wave solutions for the septic model as functions of the propagation distance. The setup and parameters are the same as in Fig. (2.26), except for M = 0.99 and F = sn. 69 Chapter 3 Application of the F-expansion technique to the Gross-Pitaevskii equation 3.1 Introduction Gross-Pitaevskii equation (GPE) is of tremendous importance in Bose-Einstein con- densation (BEC), where it describes the behavior of the condensate wavefunction (160; 161). It has been introduced independently by Gross (162; 163) and Pitaevskii (164; 165) for an unrelated problem, but has since been found of great use in BEC. In addition, it has been used in the studies of superfluidity in liquid He II, as well as pulse propagation in nonlinear (NL) optics (4). The GPE emerges in the study of the BEC from analyzing the interaction of N bosons in a pseudo-potential, which is applicable in the dilute limit, where the average spacing between the particles is much greater than the scattering length (166). The interaction between the two particles is represented by a Dirac delta function ultimately leading to a third-order nonlinearity in the GP equation. The external potential is typically taken to be quadratic in all spatial directions. Solutions to GPE are of great interest because they can be applied to a diverse array of quantum systems. Among other, solitary wave solutions (151; 167; 4) have 70 been discovered in GPE. Various, mostly numerical (168; 169), solutions to GPE have been found, prominently including localized wave solutions (170; 171). How- ever, just as in the case of the ordinary nonlinear Schro¨dinger equation (NLSE), stable exact soliton solutions to GPE exist only in (1+1)-dimensions ((1+1)D) (172; 173) and there are no known exact stable solitons in higher dimensions. In a variational and numerical treatment, Adhikari has shown that the 3D spatiotempo- ral (ST) optical solitons can be stabilized by a rapidly oscillating scattering length or the dispersion coefficient in a Kerr medium with cubic nonlinearity (174; 175). Here we present analytical traveling wave and soliton solutions to the GPE in (3+1)D, i.e. in 3 spatial dimensions and time, that were found applying the F- expansion technoique (FET). The work in this chapter was done in (145), (146), (147) and (148) . The solutions we find depict the way in which an initial traveling wave packet obeying GPE changes in time. Such solutions are necessarily transient in nature. This is the consequence of not only the equation being of the time- dependent Schro¨dinger type, but also of the fact that the coefficients in the equation are time-dependent, which is typical of BEC. Hence the solutions might diminish in time, or blow up, or oscillate, or converge to a specific spatial form. We are most interested in those solutions which oscillate or converge. The stability of these solutions will be adressed in Chapter 6. 3.2 Application of the F-expansion technique We consider GPE in (3+1)D with distributed coefficients (151): i∂tu+ β(t) 2 ∆u+ χ(t)|u|2u+ α(t)r2u = iγ(t)u. (3.1) Here t is time, ∆ = ∂2x+∂ 2 y+∂ 2 z is the 3D Laplacian, r = √ x2 + y2 + z2 is the position coordinate, and α(t) stands for the strength of the quadratic potential as a function of time. It is strictly assumed that α(t) 6= 0, otherwise the equation is reduced to the generalized nonlinear Schro¨dinger equation (NLSE), which has already been discussed in Secs. 2.2-2.3.4 and in (142; 156). The functions β, χ, and γ stand for the diffraction, nonlinearity, and gain coefficients, respectively. All coordinates in Eq. (3.1) are made dimensionless by the choice of coefficients. The key difference 71 between this system and that of 2.3.1 apart from the quadratic potential is that the longitudinal variable is now t instead of z, and the three spatial variables are now transverse. We apply the same form to the solution u as in Eq. (2.31), taking into account the change of longitudinal variables (i.e. putting the longitudinal variable as the first argument): u(t, x, y, z) = A(t, x, y, z) exp (iB(t, x, y, z)). (3.2) Substituting u into Eq. (3.1), two coupled equations for A and B are obtained: ∂tA+ 1 2 β(2∂xA∂xB + 2∂yA∂yB + 2∂zA∂zB+ A∆B) = γA, (3.3) −A∂tB + 1 2 β ( ∆A− A((∂xB)2 + (∂yB)2 + (∂zB)2) ) + χA3 + αr2A = 0. (3.4) The key change from Eqs. (2.35)-(2.36) is the presence of the extra term in Eq. (3.4) coming from the quadratic potential. To these equations we apply the balance principle and the F-expansion tech- nique, as described in 2.2. We seek the traveling wave solutions to Eqs. (3.3)-(3.4). We assume the same form for the amplitude A and the phase B as in Eqs. (2.32)- (2.34), again noting the change of the longitudinal variable: A = f1(t)F (θ) + f2(t)F −1(θ), (3.5) θ = k(t)x+ l(t)y +m(t)z + ω(t), (3.6) B = a(t)r2 + b(t)(x+ y + z) + e(t), (3.7) where f , g, k, l, m, ω, a, b, e are parameter functions to be determined, and F is a Jacobi elliptic function (JEF). Using the FET as described in 2.2, we obtain the 72 following system of algebraic or first-order ordinary differential equations: dfj dt + 3aβfj − γfj = 0, (3.8) dk dt + 2kaβ = 0, (3.9) dl dt + 2laβ = 0, (3.10) dm dt + 2maβ = 0, (3.11) da dt + 2βa2 − α = 0, (3.12) db dt + 2βab = 0, (3.13) dω dt + β(k + l +m)b = 0, (3.14) de dt + β 2 (3b2 − (k2 + l2 +m2)c2)− 3χf1f2 = 0, (3.15) f1(β(k 2 + l2 +m2)c4 + χf 2 1 ) = 0, (3.16) f2(β(k 2 + l2 +m2)c0 + χf 2 2 ) = 0, (3.17) where j = 1, 2. The constants c0, c2, c4 in Eqs. (3.8)-(3.17) are related to the elliptic modulus M of JEFs (see Table (1.1)). By solving Eqs. (3.8)-(3.17) self-consistently, one obtains a set of conditions on the coefficients and parameters, necessary for Eq. (3.1) to have exact traveling wave and soliton solutions (142). The existence of the coefficient α(t) makes the solution of Eq. (3.12) for the chirp significantly more difficult than that of the corresponding equation in the case of the generalized NLSE (142; 156). Since the solution of a(t) also determines the solutions of all other parameters, the solutions obtained here are markedly different from those in (142) and (156). Indeed, Eq. (3.12) is of the Riccati equation type, which has no analytical solutions for the general functions α(t) and β(t), although its numerical solution entails little difficulty. However, for certain choices of α(t) and β(t) it is possible to obtain exact solutions. 3.3 Solutions for proportional α and β In this section we will focus on the solutions obtained when the functions α and β are proportional. In other words, we assume α(t) = α0g(t) and β(t) = β0g(t), where 73 α0 and β0 are initial values of α and β respectively, i.e. constants, and g(t) is some function in time. In making these assumptions the Riccati equation given in (3.12) becomes separable, and therefore solvable. The final solutions to Eqs. (3.8)-(3.17) are: f1(t) = f0p 3/2 exp ( ∫ t 0 γdt), (3.18) f2(t) = ǫ √ c0 c4 f1, (3.19) k(t) = pk0, (3.20) l(t) = pl0, (3.21) m(t) = pm0, (3.22) ω(t) = ω0 − (k0 + l0 +m0)b0q, (3.23) b(t) = pb0, (3.24) e(t) = e0 + 1 2 ((k20 + l 2 0 +m 2 0)(c2 − 6ǫ √ c0c4)− 3b20)q, (3.25) where p(t) and q(t) are parameters to be determined, as well as function a(t). One can then obtain the final solution: u = f0p 3/2 exp (∫ t 0 γdt )( F (θ) + ǫ √ c0 c4 1 F (θ) ) · (3.26) exp i(a(x2 + y2 + z2) + b(x+ y + z) + e), where: θ = ω0 + (k0x+ l0y +m0z)p− (k0 + l0 +m0)b0q. (3.27) The parameter ǫ can take the values 0,±1. As a consequence of Eqs. (3.16)-(3.17), our solution method imposes an integrability condition on the coefficients: χ(t) = −β(t)c4(k20 + l20 +m20)f−20 p−1 exp ( −2 ∫ t 0 γ(t)dt ) . (3.28) Equations (3.18)-(3.28) hold for arbitrary forms of α and β, although in general a, p and q cannot be found. However, for α(t) = α0g(t) and β(t) = β0g(t) we obtain: p(t) = √ α0 α0 − 2a20β0 sech (τ(t)) , (3.29) q(t) = √ α0β0√ 2(α0 − 2a20β0) tanh (τ(t))− a0β0 α0 − 2a20β0 , (3.30) a(t) = √ α0 2β0 tanh (τ(t)), (3.31) 74 where: τ(t) = arctanh ( a0 √ 2β0 α0 ) + √ 2α0β0 t∫ 0 g(t)dt. (3.32) In (146) there was a minor error in the corresponding formula for (3.32), in that instead of g(t), β(t) was erroneously written. However, since we always used β0 = 1, this produced no additional errors. The form of the auxiliary function τ(t) naturally depends on the form of g. The parameter functions p and q place the following restriction on the solutions: α0 > 2a 2 0β0, which for positive α0 and β0 implies |a0| <√ α0 2β0 , though alternate formulas can be found for |a0| > √ α0 2β0 . We will now study the solutions obtained for two characteristic cases: constant α and β, i.e. g(t) = 1, and sinusoidal α and β, i.e. g(t) = cos (Ωt) and g(t) = sin (Ωt). In the final case α0 and β0 obviously refer to the amplitude of the functions, and not the initial values. 3.3.1 Solutions for constant α and β We now consider the special case when α and β are constants, i.e. g(t) = 1. Two distinct situations arise: either α and β are of the same sign, or they are of the opposite sign. 3.3.1.1 Case 1: α and β of the same sign We will first cover the case when α and β are of the same sign. The work in this section was done in (145). By plugging in g = 1 into Eqs. (3.29)-(3.31), we obtain the following forms for a, p and q: p(t) = est/2(1 + C) 1 + Cest , (3.33) q(t) = (1 + C)(est − 1) s(1 + Cest) , (3.34) a(t) = √ α 2β Cest − 1 Cest + 1 , (3.35) where: C = √ α 2β + a0√ α 2β − a0 (3.36) 75 and s = 2 √ 2αβ. When α and β are of the same sign, both C and s are real. The subscript 0 denotes the value of the given function at t = 0. Equations (3.34)- (3.35) hold for a0 6= √ α 2β . When a0 = √ α 2β one obtains the appropriate solution expressions by taking the limit C →∞. The final solution for u becomes: u = f0 ( est/2(1 + C) 1 + Cest )3/2 exp (∫ t 0 γdt )( F (θ) + ǫ √ c0 c4 1 F (θ) ) (3.37) · exp (i(ar2 + b(x+ y + z) + e)), where θ = ω0 + kx+ ly +mz − β(k0 + l0 +m0)b0 (1 + C)(e st − 1) s(1 + Cest) . (3.38) There are few key differences between the solutions obtained here and the ones obtained in (142). Most notably, there is no meaningful distinction in the sense of chirp vs. no chirp, i.e. between the solutions with a0 6= 0 and a0 = 0. The value of a0 = 0 does not entail any special status. Instead, it is the value of √ α 2β that is of some importance. For a0 > − √ α 2β , a converges to √ α 2β as t increases, while for a0 = ± √ α 2β it stays constant and for a0 < − √ α 2β there are singularities in the parameter a. For a0 > − √ α 2β the functions k, l, m and b all converge to 0 and ω and e converge to constant values that depend on a0. Figure 3.1: Decaying bent solitary wave solutions to GPE as functions of time, for b0 = 1. Intensity |u|2 for: (a) F = tanh and (b) F = sech presented as functions of k0x + l0y +m0z and t. Other parameters are: β = 1, α = 1, γ(t) = −0.05, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. 76 Figure 3.2: Decaying straight soliton solutions to GPE as functions of time. The setup and parameters are the same as in Fig. (3.1) except for b0 = 0. Figure 3.3: Decaying traveling wave solutions, given in terms of JEFs for: (a) F = sn and (b) F = cn. The setup and parameters are the same as in Fig. (3.2), except for M = 0.99. 77 By inspecting Eq. (3.42) one can see that the long-time behavior of the general solution crucially depends on the coefficient γ(t). Although γ(t) is described as the linear gain or loss in the system, the value of γ = 0 does not exert any special bearing on the solutions, similar to the value of a0 = 0 for chirp. Figures (3.1)-(3.3) depict decaying solutions for a small negative value of γ. The solutions for γ = 0 are also decaying and are practically identical to those in Fig. (3.1). The critical value of γ for the appearance of solitons or waves as t increases is γ = 3p/4. Thus, if γ is constant, only for γ = 3p/4 can one see stable solitons or waves evolving as t → ∞. If γ > 3p/4, the solutions blow up, while if γ < 3p/4, the solutions diminish. Hence, to observe stable solitons asymptotically, one needs gain in the system. To see periodically changing, i.e. breathing, solitons in the case of constant α and β, one needs γ in the form of γ(t) = 3p/4 + γ1(t), where γ1(t) is some periodic sign-changing real function. Figure 3.4: Bent soliton solutions as functions of time for: (a) F = sn and (b) F = cn. The setup and parameters are the same as in Fig. (3.1), except for γ(t) = 3/ √ 2, the critical value of γ. The caveat to the analysis just presented is contained in Eq. (3.28), which explicitly connects χ(t) with γ(t). Thus, the long-time behavior of the nonlinearity coefficient χ is also tied to γ(t). For constant γ, the critical value now is γ = p/4. If γ < p/4, χ(t) diminishes with time, if γ > p/4, it blows up. Taken together with the result of the previous paragraph, it appears that the most interesting interval of γ for the long-time behavior is p/4 ≤ γ(t) ≤ 3p/4. There, the solutions decrease in time, while the nonlinearity coefficient increases. This statement reflects the 78 Figure 3.5: Straight soliton solutions as functions of time. The parameters are the same as in Fig. (3.4) except for b0 = 0. Figure 3.6: Traveling wave solutions in terms of JEFs. The parameters are the same as in Fig. (3.5), except for M = 0.99. (a) F = sn and (b) F = cn. 79 difficulties in obtaining stable solitons in the multidimensional GP equation with constant coefficients. It is another reflection of the known difficulties with the wave stability and collapse in multidimensional NL Schro¨dinger equation (48). Hence, to observe long-lived solitons, a delicate engineering in the form of γ(t) is necessary. As in the case of the NLSE, the solutions introduced by Eqs. (3.26)-(3.7) describe spatially extended spatio-temporal solitons for M = 1. 3.3.1.2 Case 2: α and β of opposite sign The case where α and β are of opposite sign is of much interest to the scientists studying the GP equation because in this case the solutions simulate an attractive quadratic potential, which is of a much larger relevance for experiments done on Bose-Einstein condensates, whereas when α and β of the same sign, the solutions simulate the effect of a repulsive quadratic potential. We define α′ = −α. After solving Eqs. (3.8)-(3.17), one obtains the following solutions for a, p and q: p(t) = 1 cos (√ 2α′βt ) + a0 √ 2β α′ sin (√ 2α′βt ) , (3.39) q(t) = √ α′ 2β a0 √ 2α′β + α′ cot (√ 2α′βt ) , (3.40) a(t) = √ α′ 2β  a0 − √ α′ 2β tan (√ 2α′βt ) √ α′ 2β + a0 tan (√ 2α′β )   . (3.41) The final solution for u is then: u = f0   √ α′ 2β a0 √ 2α′β + α′ cot (√ 2α′βt )   3/2 exp (∫ t 0 γdt ) · ( F (θ) + ǫ √ c0 c4 1 F (θ) ) exp (i(ar2 + b(x+ y + z) + e)), (3.42) where: θ = ω0 + kx+ ly +mz − β(k0 + l0 +m0)b0 √ α′ 2β a0 √ 2α′β + α′ cot (√ 2α′βt ) . (3.43) Unlike the solutions in Case 1, which dissipate in time, the solutions in Case 2 tend to converge and blow up. In this case, however, no amount of loss will 80 rescue the solutions from collapsing as parameter p, on which the amplitude of the solutions depends, inevitably reaches 0. For a0 = 0, The solutions blow up at time T/4, where T = 2π/( √ 2α′β). To avoid singularities we show all graphs to almost, but not including, the time at which the singularity is achieved. In Fig. (3.7) we see see the form of the solutions for both the bright and dark solitary waves and the effect of parameter b0 on the bright solitary waves. The parameter b0 ultimately determines the direction of the solitary wave. Unlike the case of the NLSE where the wave is uniformly travelling in one direction for constant β, here the wave travels sinusoidally, but blows up after only a quarter of a period. Figure 3.7: Solitary wave solutions to GPE as functions of time. Values of param- eters are: β = 1, α = −1, γ(t) = −0.05, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. Intensity |u|2 for (a) F = tanh, b0 = 2 and for the remaining graphs F = sech with: (b) b0 = 2, (c) b0 = 0 and (d) b0 = −2 presented as functions of k0x+ l0y +m0z and t, where t is shown from 0 to T/4− 0.01. The effects of chirp are similar to Case 1. There is no qualitative change in the 81 form of the solutions, and the presence of a0 does not prevent the solution from blowing up. As a0 increases and the sin term in the denominator of p becomes more prominent, the time at which the solutions blow up converges to T/2. In addition, the presence of a0 causes a decrease in amplitude with respect to the starting point. We see these effects in Fig. (3.8). Unlike in Case 1, or indeed in the case of the NLSE described in Sec. 2.3.1, the initial value of the chirp a0 can be made as large as desired without running into (new) singularities. Figure 3.8: Solitary wave solutions to the GPE as functions of k0x+ l0y +m0z and t, where t is shown from 0 to T/2 − 0.01. The setup and parameters are the same as in Fig. (3.7)(d) except for: (a) a0 = 0.1, (b) a0 = 0.5, (c) a0 = 1 and (d) a0 = 2. In Fig. (3.9) we see the travelling wave solutions for Case 2. As in Case 1, the solutions are not periodic with respect to the transverse variable. Unlike Case 2, the travelling waves converge to a single point, whose location depends on b0. 82 Figure 3.9: Traveling wave solutions to the GPE as functions of k0x+ l0y+m0z and t, where t is shown from 0 to T/2 − 0.01. The setup and parameters are the same as in Fig. (3.7) except for M = 0.99, a0 = 1 and: (a) F = sn, b0 = 0, (b) F = cn, b0 = 0, (c) F = sn, b0 = 10 and (d) F = cn, b0 = 10. 83 3.3.2 Solutions for sinusoidal α and β We now present and analyze the solutions to GPE when both the diffraction and the quadratic potential are sinusoidal functions of time. As in Sec. 3.3.1, there will be two distinct cases: either α0 and β0 are of the same sign, in which case the two functions α(t) and β(t) are in phase, or α0 and β0 are of the opposite sign, in which case the two functions α(t) and β(t) have a phase difference of 180◦. We attempted to solve other cases, but unsuccessfully. If arbitrary phase shifts are added to α and β, i.e. α(t) = α0 sin(Ωt+φ1) and β(t) = β0 sin(Ωt+φ2), very different solutions are found, often without a closed analytical form. Indeed, Equation (3.12) is separable only for φ1 = 0 and φ2 = π. If one parameter is a sine function and the other is a constant, the solutions collapse rapidly. The GP equation will be solved in each case Eq. (3.1) for two sub-cases: α(t) = α0 sin(Ωt), β(t) = β0 sin(Ωt) (g(t) = sin(Ωt)) and α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt) (g(t) = cos(Ωt)). 3.3.2.1 Case 1: α0 and β0 of the same sign We first consider that α0 and β0 are of the same sign. This case was analyzed for the first time in (146). The form of parameter functions p and q, as well as the solution for a, differ in the two sub-cases. For the first case, α(t) = α0 sin(Ωt), β(t) = β0 sin(Ωt), we have: a = √ α0 2β0 tanh ( arctanh ( a0 √ 2β0 α0 ) + √ 2α0β0 (1− cos (Ωt)) Ω ) , (3.44) p = sech ( arctanh ( a0 √ 2β0 α0 ) + √ 2α0β0 (1− cos (Ωt)) Ω ) · (3.45) √ α0 α0 − 2a20β0 , q = tanh ( arctanh ( a0 √ 2β0 α0 ) + √ 2α0β0 (1− cos (Ωt)) Ω ) · (3.46) √ α0β0√ 2(α0 − 2a20β0) − a0β0 α0 − 2a20β0 . 84 For the second case, α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt), we have: a = √ α0 2β0 tanh ( arctanh ( a0 √ 2β0 α0 ) + √ 2α0β0 sin (Ωt) Ω ) , (3.47) p = √ α0 α0 − 2a20β0 sech ( arctanh ( a0 √ 2β0 α0 ) + √ 2α0β0 sin (Ωt) Ω ) , (3.48) q = √ α0β0√ 2(α0 − 2a20β0) tanh ( arctanh ( a0 √ 2β0 α0 ) + √ 2α0β0 sin (Ωt) Ω ) (3.49) − a0β0 α0 − 2a20β0 . We present in this section a few typical examples of solutions for both sub-cases, α(t) = α0 sin(Ωt), β(t) = β0 sin(Ωt) and α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt). The initial conditions (sin 0 = 0, cos 0 = 1) produce a crucial difference in the chirp parameter a, which in turn affects both p and q. Another important point to note is that while the traveling wave solutions are periodic in time, they are not periodic along the transverse variable k0x+ l0y +m0z, in contrast to the solutions found in Sec. 2.3.1. Like in Sec. 3.3.1, the initial value of the chirp is not of much importance, i.e the solutions remain qualitatively the same. The only major change is a shift of all parameters a, p and q, which causes a shift in the graphs along the transverse variable, and a decrease in the magnitude (for positive a0) which causes a narrowing of the peaks. Figure (3.10) presents the sine case and Fig. (3.10) the cosine sub-case of the solitary wave. Figure (3.12) presents the traveling wave solutions for both the sine and the cosine sub-case. For a better perspective, a small loss (γ = −0.05) is included in all the figures. Without it, the waves still breathe, the amplitude remains constant. Note the influence of the parameter b0, which causes the solitons to wiggle. Despite the apparent complexity of these solutions they yield relatively straight- forward and elegant spatiotemporal breathing soliton solutions for both F = sn and F = cn, i.e. both dark and bright solitons, without any need for γ. In other words, unlike the solutions in Sec. 3.3.1 which required a positive value of gain to form soli- tary waves, the signals here stay at the same peak intensity for γ = 0, but breathe. This is apparent from the fact that p and q, as well as a, the parameter functions on which all other variables in the final solution depend, are periodic functions of 85 Figure 3.10: Soliton solutions to Gross-Pitaevskii equation as functions of time, for the sine case: α(t) = α0 sin(Ωt), β(t) = β0 sin(Ωt). Intensity |u|2 presented as functions of k0x+ l0y +m0z and t for: (a) F = tanh, b0 = 0, (b) F = sech, b0 = 0, (c) F = tanh, b0 = 1 and (d) F = sech b0 = 1. Other parameters are: α0 = 1, β0 = 1, Ω = 1, γ(t) = −0.05, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0. 86 Figure 3.11: Soliton solutions to Gross-Pitaevskii equation as functions of time for the cosine case: α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt). Other parameters are the same as in Fig. (3.10) 87 Figure 3.12: Traveling wave solutions to Gross-Pitaevskii solutions in terms of JEFs for the sine and the cosine case. The parameters for (a) and (b) are the same as in Fig. (3.10) (a) and (b), and the parameters in (c) and (d) are the same as in Fig. (3.11) (a) and (b), except for M = 0.99. For figures (a) and (c) we have F = sn and for figures (b) and (d) we have F = cn. 88 time in both cases. Note also that the width of the solitary solutions in Sec. 3.3.1 was increasing in time, because of the positive value of γ. 3.3.2.2 Case 2: α0 and β0 of the opposite sign In Case 2, we cover solutions in which α0 and β0 of the opposite sign. The formulas for f1, f2, b, k, l, m, ω and e are the same as in Eqs. (3.18)-(3.25). The forms for a and parameters p and q for the case α(t) = α0 sin(Ωt), β(t) = β0 sin(Ωt), are: a = √ α0 2β′0 tan ( arctan ( a0 √ 2β′0 α0 ) − √ 2α0β′0 (1− cos (Ωt)) Ω ) , (3.50) p = √ α0 α0 + 2a20β ′ 0 sec ( arctan ( a0 √ 2β′0 α0 ) − √ 2α0β′0 (1− cos (Ωt)) Ω ) ,(3.51) q = tan ( arctan ( a0 √ 2β′0 α0 ) − √ 2α0β′0 (1− cos (Ωt)) Ω ) · (3.52) √ α0β′0√ 2(α0 + 2a20β ′ 0) + a0β ′ 0 α0 + 2a20β ′ 0 , where β′0 = −β0. For the second case, α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt), we have: a = √ α0 2β′0 tan ( arctan ( a0 √ 2β′0 α0 ) − √ 2α0β′0 sin (Ωt) Ω ) , (3.53) p = √ α0 α0 + 2a20β ′ 0 sec ( arctan ( a0 √ 2β′0 α0 ) − √ 2α0β′0 sin (Ωt) Ω ) , (3.54) q = √ α0β′0√ 2(α0 + 2a20β ′ 0) tan ( arctan ( a0 √ 2β′0 α0 ) − √ 2α0β′0 sin (Ωt) Ω ) (3.55) + a0β ′ 0 α0 + 2a20β ′ 0 . Here we present the results for just the case α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt). Apart from the shape of the peaks and their ultimately large value for the same set of parameters, as seen in Fig. (3.13), the results for Case 2 are very similar to those in Case 1. The results for the α and β being sine functions are also similar, though singularities occur for different sets of values of α0 and β ′ 0. For example, for the same set of parameters as those in Fig. (3.13) singularity occurs in the sine case. 89 Figure 3.13: Soliton solutions to Gross-Pitaevskii equation as functions of time for the cosine case: α(t) = α0 cos(Ωt), β(t) = β0 cos(Ωt). Other parameters for figures (a) and (b) are the same as in figures Fig. (3.10)(c) and (d), respectively, except for b′0 = −b0 = 1. 3.4 Other systems with solvable Ricatti equations We now turn to solving other systems for which one can solve the corresponding Ricatti equation (3.12). The solution procedure employed here is different from the procedure in (145) and (146). The procedure there requires that α(t) and β(t) be of the same signs, whereas here they have to be of the opposite sign. Also, the notation introduced in Sec. 3.4.1, for simplicity’s sake, differs somewhat from notation elsewhere in the Thesis, and thus the two should not be confused. 3.4.1 Solution method The Riccati equation (RE) has the following form: y′ = P (x) +Q(x)y +R(x)y2, (3.56) which can be considered as the lowest order nonlinear approximation to the deriva- tive of a function in terms of the function itself (176). It is assumed that y, P , Q and R are real functions of the real argument x. It is well known that the general so- lution to the Riccati equation is not available, and only special cases can be treated 90 (177; 178; 179; 180; 181; 182). Even though the equation is nonlinear, similar to the second order inhomogeneous linear ODEs, one only needs a particular solution to find the general solution. In a standard manner Riccati equation can be reduced to a second-order linear ODE (177) or to a Schro¨dinger equation (SE) of quantum mechanics (183). In fact, Riccati equation naturally arises in many fields of quantum mechanics, in particular, in quantum chemistry (184), in the use of the Wentzel-Kramers-Brillouin (WKB) approximation (185) and in SUSY theories (187). Recently, methods for solving the Gross-Pitaevskii equation (GPE) arising in Bose-Einstein condensates (BECs) (172; 173) based on the Riccati equation were introduced. Our objective is to find new solutions of Riccati equation by utilizing relations between the coefficient functions P (x), Q(x) and R(x) for which the above equation can be solved in closed form. It is well known that any equation of the Riccati type can always be reduced to the second order linear ODE: u′′ − ( Q(x) + R′(x) R(x) ) u′ + P (x)R(x)u = 0 (3.57) by a substitution: y = − u ′ uR . (3.58) It is also known that if one can find a particular solution yp to the original equation, then the general solution can, in a well-known procedure, be written as (147): y = yp + 1 w , (3.59) where w is the general solution of an associated linear ODE: w′ + (Q(x) + 2R(x)yp)w +R(x) = 0 (3.60) which does not contain P (x). Solving this equation we get (147): w = w0e −φ(x) − e−φ(x) ∫ x x0 R(ξ)eφ(ξ) dξ, (3.61) where: φ(x) = ∫ x x0 (Q(ξ) + 2R(ξ)yp) dξ. (3.62) 91 It is clearly seen from the relation above that w0 = 1 y0−yp0 . The general solution is therefore given by (147): y = yp + e φ(x) ( 1 y0 − yp0 − ∫ x x0 R(ξ)eφ(ξ) dξ )−1 . (3.63) Equation (3.56) cannot be solved in closed form for arbitrary functions P (x), Q(x) and R(x). However, if certain relations exist between these functions, then the above equation can be transformed into a second order linear ODE, which can be easily solved in certain special cases, which we will cover. For the sake of making our calculations clearer, we make the following two substitutions: a(x) = − ( Q+ R′ R ) , (3.64) and: b(x) = P (x)R(x). (3.65) Parameters a(x) and b(x), which we will use only in this section, should not be con- fused with the chirp function a(t) and with the function b(t), both used extensively elsewhere in the Thesis. Equation (3.57) can thus be re-written as: d2u dx2 + a(x) du dx + b(x)u = 0. (3.66) We consider an arbitrary function of x, z ≡ f(x), which we choose to be a new independent variable. The substitution looks arbitrary, but it will be made more specific later in the text. We compute the first and second derivatives of u with respect to x, but now in terms of the new independent variable z: du dx = du dz dz dx , (3.67) d2u dx2 = d2z dx2 du dz + ( dz dx )2 d2u dz2 . (3.68) We plug the last results into the Equation (3.66), to get:( dz dx )2 d2u dz2 + ( d2z dx2 + a(x) dz dx ) du dz + b(x)u = 0. (3.69) Finally, dividing by ( dz dx )2 , we obtain (147): d2u dz2 + ( d2z dx2 + a(x) dz dx( dz dx )2 ) du dz + ( b(x)( dz dx )2 ) u = 0 (3.70) 92 provided dz/dx is not equal to 0. Hence, we define parameters A and B as follows: A = d2z dx2 + a(x) dz dx 2 ( dz dx )2 , (3.71) B = b(x)( dz dx )2 , (3.72) to obtain: d2u dz2 + 2A du dz +Bu = 0, (3.73) The obtained equation can easily be solved in closed form if A and B are either constants (147) or if they are some special functions for which the closed-form so- lutions to (3.73) are known. In this Thesis we consider only the two special cases, namely when A and B are positive constants, or when A = 0 and B is an arbitrary function B(x). The constants A and B should not be confused with the amplitude and phase of u given in Eq. (2.5) and elsewhere. If b(x) is positive, by considering the coefficient of u we define z to be the following function: z = z0 + s ∫ x x0 √ b(ξ) B dξ, (3.74) where s = ±1. The requirement that b(x) is positive is equivalent to the condition that the product P (x)R(x) is positive. To simplify bookkeeping, let c = b/B. In that case, we have the following relations: dz dx = sc1/2, (3.75) d2z dx2 = c′ 2sc−1/2 . (3.76) From (3.75) it is clear that dz/dx cannot be equal to 0. Now we compare the coefficients of du/dz and use relations (3.75) and (3.76) to get: c′ 2sc1/2 + asc1/2 − 2Ac = 0, (3.77) or ( b B )′ + 2a ( b B ) − 4As ( b B )3/2 = 0. (3.78) At this point it is more convenient to consider the two cases separately. 93 3.4.1.1 Case 1: A and B are constants If (3.73) has constant coefficients 2A and B, then it is easily solvable in closed form. This means: b′ + 2ab− 4sA√ B b3/2 = 0, (3.79) in other words: b′(x) + 2a(x)b(x) (b(x))3/2 = 4sA√ B . (3.80) Substituting back the original expressions for a(x) and b(x), we get the final result: (P (x)R(x))′ − 2 (Q(x) +R′(x)/R(x))P (x)R(x) (P (x)R(x))3/2 = 4sA√ B . (3.81) Note that only if the condition (3.81) is satisfied, can the general solution be found using this method. On the other hand, when the condition is not satisfied, this does not mean that the general solution cannot be found. In fact, most of the special cases of Riccati equation with known solutions do not satisfy the relation obtained (182). Also note that even though Eq. (3.79) is nonlinear, it can readily be solved. Now we proceed to solve Eq. (3.73). The general solution is given by: u(x) = c1e λ1z + c2e λ2z, (3.82) where z is the function defined in (3.74), c1 and c2 are some initial values, and λ1 and λ2 are the roots of the characteristic polynomial λ 2 + 2Aλ+B = 0, given by: λ1,2 = −A± √ A2 − B. (3.83) We assume first that A2 ≥ B > 0, so that both lambdas are real and negative. We need only a particular solution of (3.73), so we consider only up = e λz, where λ is any of the roots to the polynomial. From the substitution done in (3.57), namely y = −u′/(uR(x)), we find the particular solution to be: yp = − sλ√ B √ P (x) R(x) . (3.84) Finally, we plug yp into the expression for the general solution of Riccati equation, to find: y = − sλ√ B √ P (x) R(x) + eφ(x)   1 y0 + sλ√ B √ P (x0) R(x0) − ∫ x x0 R(ξ)eφ(ξ) dξ   −1 . (3.85) 94 Note that we have substituted yp0 by its value. To recapitulate, here A and B are two arbitrary constants satisfying A2 ≥ B > 0, λ is one of the roots of the characteristic polynomial, y0 is the value for y at x0, and: φ(x) = ∫ x x0 ( Q(ξ)− 2sλ√ B √ P (ξ)R(ξ) ) dξ (3.86) is the integrating exponent. In the case of A2 < B we obtain trigonometric solutions with damping. The general solution for u is: u(x) = e−Az(c1 cosωz + c2 sinωz), (3.87) where ω = √ B − A2. Taking c1 = 1 and c2 = 0, since we only need a particular solution, we obtain: yp = sA√ B (1− ω tanωz) √ P (x) R(x) . (3.88) Note that an explicit dependence of yp on z occurs, unlike in the case A 2 ≥ B > 0. Bearing in mind the dependence of z on x given in Eq. (3.74) and the dependence of φ on yp given in Eq. (3.62), the final formula for y for the case A 2 < B becomes prohibitively complex, involving a double integral where the inner one is embedded inside a trigonometric function. In this Thesis we will only cover the results for A2 ≥ B > 0. Below we apply the general results of this section to some specific examples. 3.4.1.2 Case 2: A = 0 and B = B(x) When A = 0, (3.77) reduces to the simple equation c′ = −2ac. Solving for c, and taking into account that c = b/B, we get the simple relation: b B = ( b B ) 0 exp ( −2 ∫ x x0 a dx ) , (3.89) where a(x) and b(x) are defined in Eqs. (3.64)-(3.65), and B(x) is still an arbitrary function. Note that (3.73) now becomes: d2u dz2 +B(z)u = 0, (3.90) where z is given by (3.74). When B(z) is chosen as B(z) = B0+B1(z), then the last equation becomes equivalent to the Schro¨dinger equation of quantum mechanics: ψ′′ + 2m ~2 (E − V )ψ = 0, (3.91) 95 describing the wave function ψ = u(z) of a particle of mass m moving in a potential: V = −~ 2B1(z) 2m (3.92) with an energy eigenvalue: E = ~ 2B0 2m (3.93) (186). There are many specific potentials V for which the solutions ψn and the energies En in the above equation are known, where n denotes some set of quantum numbers. Therefore, one can choose B(z) such that the solutions un(z) can be found. If un(z) are known, and hence un(x), then the solutions yn to Riccati equation can be easily written down from the substitution yn(x) = − u ′ n unR (3.94) mentioned above. This in fact gives rise to various solutions of the various special cases of Riccati equation. 3.4.2 Application of the solution method We now apply the results of the previous section to the GPE, as given in Eq. (3.1). The solution to the system is given by Eqs. (3.18)-(3.25). The key challenge is finding the solution for the chirp function a(t) as well as parameters p and q, although for b0 = 0 the parameter q is not needed for plotting |u|2. Again, the functions a(t) and b(t) given in Eq. (3.31) and Eq. (3.24) respectively should not be confused with the functions a(x) and b(x) defined in Eqs. (3.64)-(3.65). The equation for a(t), (3.31), will be restated here: da dt + 2β(t)a2 − α(t) = 0. (3.95) To this equation we apply the method developed in Sec. 3.4.1. We take parameters A and B to be constant, i.e. we only cover Case 1 in this Thesis. Put in the form of the original Riccati equation, the coefficients are: P (t) = α(t), (3.96) Q(t) = 0, (3.97) R(t) = −2β(t). (3.98) 96 We write down relation (3.81) between α and β for which (3.95) is solvable in closed form: αβ′ − α′β (−αβ)3/2 = 4 √ 2sA√ B . (3.99) The prime is now the derivative with respect to t. Equation (3.99) can be manipu- lated to become a simple differential equation for −α/β:(− α β )′ (− α β )3/2 β = 4 √ 2sA√ B . (3.100) Solving this equation, one finds:√ −β α = √ −β0 α0 − 2 √ 2sA√ B ∫ t 0 β dt. (3.101) Now one can write down the solution for a(t) from (3.85), provided the above con- dition is satisfied: a(t) = − sλ√ B √ − α(t) 2β(t) + eφ(t)   1 a0 + sλ√ B √ − α0 2β0 + 2 ∫ t 0 β(τ)eφ(τ) dτ   −1 , (3.102) where: φ(t) = −2 √ 2sλ ∫ t 0 √ −α(τ)β(τ) dτ/ √ B. (3.103) Note that the negative sign in the square root indicates that α and β have to be of the opposite signs, which is consistent with the requirement that the original function b(x) is positive. Hence, as long as the ratio of the diffraction coefficient to the strength of the parabolic potential can be made to satisfy (3.100), one can write down the exact solutions to GPE. It should be mentioned that these functions are the material parameters in BECs that are accessible to experimental manipulation. Our solution method for the GPE requires that β be proportional to χ, and χ in turn be proportional to the s-wave scattering length (188). To validate our proposed solution method, we present a couple of examples in which β, and hence the scattering length, are given by some representative functions of time. In all the examples we determine the corresponding chirp functions a(t), from which one can write down the exact solutions of the GPEs in question (145). To avoid singularities that are likely to appear in α(t) and a(t) we choose s to be −1. Note that the appearance of singularities is not detrimental to our method or to the theory of BECs based on GPE, because that model is known to be valid only on a limited 97 time interval. In the subsequent examples, the calculations and graphs for α, φ and a(t) were done in (147), whereas the remaining work, including the graphs for u, was done in (148). The three examples we will cover are possible to implement in real physical systems and, as far as we know, are new solutions which were not mentioned in the literature before. 3.4.2.1 Example 1: β = 1 2 (e−δt + 1) We consider first the case when β is an exponential function of time, β(t) = 1 2 (e−δt+ 1), where δ is some arbitrary parameter. This function describes a smooth change in β(t) from 1 to 1/2. First, (3.101) is solved for α, to obtain: α(t) = − 1 + e −δt 2 ( 1 + √ 2(1−e−δt+δt) δ )2 . (3.104) Then one finds φ : φ(t) = δt+ ln ∣∣∣∣∣ δ−√2 + eδt (√2 + δ +√2δt) ∣∣∣∣∣ . (3.105) Taking α0 = −1, A = B = 1, and performing the calculations, we obtain the following solution for a: a(t) = −δ 2− 2e−δt +√2δ + 2tδ − δ √ 2eδt(−√2 + eδt (√2 + δ +√2tδ)) ζ(t) , (3.106) where ζ(t) = δt+ ln ∣∣∣∣ δ−√2 + etδ(√2 + δ +√2tδ) ∣∣∣∣− 21 +√2a0 . Although this solution looks complicated, it allows simple expressions in the limit δ → 0, when β becomes constant. Figure (3.14) presents some representative cases of α and a functions for different values of δ. Finally, the parameter p is calculated by solving Eq. (3.24) for Eq. (3.106): p = − 2 5/4e− δt 2 √ −2+eδt(2+( √ 2+2t)δ) δ −2√2 +√2tδ + 2tδa0 + ln ∣∣∣∣ δ−√2+eδt(√2+δ+√2tδ) ∣∣∣∣ (√2 + 2a0) . (3.107) Figure (3.15) shows how the solution of GPE looks like for this case. One may note that after an initial rapid change, the pulse settles into a slowly evolving bright solitary solution. Similarly to solutions in Sec. 3.3.1 there is a continuous reshaping 98 Figure 3.14: Graphs of parameters for Example 1: (a) α(t), (b) a(t) for a0 = 0 and (c) a(t) for a0 = 1, for δ = 0.01, 0.1, 1, 10 (top to bottom). Figure 3.15: Intensity distribution |u|2 for the solution of Example 1. (a): no gain/loss. (b): γ = −0.05. Here F = sech. Other parameters are: a0 = f0 = k0 = l0 = m0 = ω0 = 1, b0 = ǫ = 0 and δ = 5. of the wave, which is relatively mild in this example and Example 3, but far more dramatic in Example 2. For Example 1, it can be inferred that the intensity linearly grows with time, even when there is no gain imposed on the system. To prevent the intensity from be- coming arbitrarily large, some loss should be added to the system. In Fig. (3.15)(b) this is achieved by choosing γ = −0.05. 3.4.2.2 Example 2: β = ∑N n=0 βnt n Next we consider the case when β is some power series of the form ∑N n=0 βnt n, where β0 6= 0. We go through the same procedure and solve (3.101) for α, to get: α(t) = − ∑N n=0 βnt n( 1 + 2 √ 2 ∑N n=0 βn tn+1 n+1 )2 . (3.108) 99 Then we find φ to be: φ(t) = ln 1∣∣1 + 2√2∑Nn=0 βn tn+1n+1 ∣∣ . (3.109) Again, taking α0 = −1, A = B = 1, and performing the calculations, we arrive at the following closed-form solution: a(t) = 2 √ 2a0 − (a0 √ 2 + 1) ln ∣∣1 + 2√2∑Nn=0 βn tn+1n+1 ∣∣( 1 + 2 √ 2 ∑N n=0 βn tn+1 n+1 )( 2 √ 2 + (2a0 + √ 2) ln ∣∣∣1 + 2√2∑Nn=0 βn tn+1n+1 ∣∣∣ ) . (3.110) These solutions for α and a are plotted in Fig. (3.16). Note that by choosing different parameters βn and letting N →∞ one can obtain closed-form expressions for different functions β(t). Figure (3.17) presents the case with β = cos(Ωt). Figure 3.16: Graphs of parameters for Example 2: (a) α(t), (b) a(t) for a0 = 0, and (c) a(t) for a0 = 1. Parameters: N = 0, 1, 2, 3, 4 (top to bottom at t = 0.5 for α, bottom to top at t = 3 for a), βn = 1. Figure 3.17: Graphs of parameters for β(t) = cos(Ωt): (a) α(t), (b) a(t) for a0 = 0 and (c) a(t) for a0 = 1. Here Ω = 6, 7, 8, 9, 10. Curves with higher peaks correspond to lower values of Ω. We now study the special case: β = ∞∑ n=0 (−1)n (2n)! (Ωt)2n = cos(Ωt), (3.111) 100 which was covered in (148). We find α to be: α = − cos(Ωt)( 1 + 2 √ 2 sin(Ωt) Ω )2 . (3.112) The chirp function then has the following form (147): a(t) = 2 √ 2a0 − (√ 2a0 + 1 ) ln ∣∣∣2√2 sin(Ωt)Ω + 1∣∣∣( 2 √ 2 sin(Ωt) Ω + 1 )(( 2a0 + √ 2 ) ln ∣∣∣2√2 sin(Ωt)Ω + 1∣∣∣+ 2√2) . (3.113) The expression for p becomes: p = 2 √ 2 + 4 √ 2 sin(Ωt) Ω 2 √ 2 + ln ∣∣∣1 + 2√2 sin(Ωt)Ω ∣∣∣ (√2 + 2a0) . (3.114) Next, we find q(t): q = er ( 2Ω2 (( 2a0 + √ 2 ) ln (ξ2(t)) + 2 √ 2 ) Ei (2 ln (ξ2(t))− r)− e−rξ1(t) ) Ω2 (( a0 ( 2a0 (√ 2a0 + 3 ) + 3 √ 2 ) + 1 ) ln (ξ2(t)) + 4a0 ( a0 + √ 2 ) + 2 ) + ξ3(t), (3.115) where: ξ1(t) = 8( √ 2a0 + 1)Ω sin(Ωt)− 4(2a0 + √ 2) cos(2Ωt) + (3.116) (2a0 + √ 2)(Ω2 + 4), ξ2(t) = 2 √ 2 sin(Ωt) Ω + 1, (3.117) ξ3(t) = −4e− 4√ 2a0+1Ei ( 4√ 2a0+1 ) + √ 2a0 + 1 √ 2 (√ 2a0 + 1 )2 , (3.118) r = − 4√ 2a0 + 1 , (3.119) and Ei(x) = ∫ x −∞ et t dt is the exponential integral. Note that here a proper choice of Ω had to be made, to ensure that the solutions do not blow up. Figure (3.18) displays the behavior of the solution. The obtained solutions greatly resemble those in Sec. 3.3.2. The choice of periodic β produces a breathing localized solution, i.e. a breather. The solution looks like a regular breather when the parameter function b of the solution is equal to 0, but wiggles back and forth keeping an asymmetric profile when b 6= 0. Such solutions propagate stably, with a periodic change in the profile. It can be clearly observed that the addition of the parameter b0 causes the periodic change in the soliton’s direction and shape. The amplitude of the soliton does not change when there is no gain or loss added to the system. 101 Figure 3.18: Intensity distribution in Example 2, with b0 = 0 (left), and b0 = 5 (right). Here F = sech. In both cases, Ω = 8. Other parameters: a0 = f0 = k0 = l0 = m0 = ω0 = 1, γ = ǫ = 0. 3.4.2.3 Example 3: β = β˜ ( 1− D B1t−B0 ) Finally, we consider the case when β is of the form: β = β˜ ( 1− D B1t−B0 ) . This form is dictated by the dependance of the scattering length on the magnetic field near the Feshbach resonance of cold BEC atoms (188). The magnetic field B(t) = B1t, not to be confused with the function B(x) from Sec. 3.4.1, is assumed to be linearly ramped in time near the resonance field B0. The parameter D stands for the width of the resonance. Such a dependence is relevant not only for for theoretical physics (188), but also, most importantly, experimental physics (189). The closed-form solution is again readily obtained. However, this time it includes integrals that cannot be evaluated in terms of elementary functions. The results for α, φ, and a are as follows: α(t) = −β˜ 1− D B1t−B0( 1− 2√2β˜t+ 2 √ 2β˜D B1 ln ∣∣B1t−B0 B0 ∣∣)2 , (3.120) φ(t) = 2 √ 2 β˜ ∫ t 0 1− D B1τ−B0 1− 2√2β˜τ + 2 √ 2β˜D B1 ln ∣∣B1τ−B0 B0 ∣∣ dτ, (3.121) a(t) = 1√ 2− 2β˜t+ 4β˜D B1 ln ∣∣B1t−B0 B0 ∣∣ + (3.122) eφ(t) √ 2 a0 √ 2−1 + 2β˜ ∫ t 0 ( 1− D B1τ−B0 ) eφ(τ) dτ . The parameter p unfortunately cannot be found in terms of simple elementary func- tions. As a consequence, numerical evaluation of p, in the absence of a closed 102 Figure 3.19: Intensity distribution in case 3, with no gain/loss (γ = 0). Here, F = sech. Other parameters D = 10, B0 = −1, a0 = B1 = f0 = ω0 = e0 = k0 = l0 = m0 = β˜ = 1, and b0 = 0. formula, had to be used in our calculations. We will not deal with the singularities that result from the resonance form of β. Therefore, we will choose parameters B1 and B0 such that the denominator remains finite. This can be accomplished, for example, if we choose B1 = 1, and B0 = −1. We choose to keep the closed-form solutions in integral form, and visualize these solutions instead. The behavior of the solution is shown in Fig. (3.19). Notice that the parameter values are properly chosen such that the solution does not blow up. The solution starts from small initial values, but rapidly grows and then continuously attenuates, the reason being the presence of resonance in the diffraction coefficient. For other values of parameters the solution might collapse. If the parameter b0 were chosen different from 0, the same solitary wave would curve and change direction. However, the wave would not wiggle back and forth in this case, as in Example 2, because the diffraction function is not periodic in time. 103 Chapter 4 Application of the F-expansion technique to the NLSE in a linear electric field 4.1 Introduction The NLSE with a cubic nonlinearity and linear potential has also been of great interest to various fields of physics. The topic was first covered in (190), where detailed solutions were offered for the (1+1)D case for constant diffraction and constant linear potential in time. Then in (114) and (115) Hirota method was used to find both one- and two-soliton solutions in (1+1)D. Finally, in (191), (1+1)D solutions without chirp for time-dependent linear potentials are developed using the F-expansion technique. Here we present analytical traveling wave and soliton solutions to the NLSE in (3+1)D with cubic nonlinearity and time-dependent linear external potential. The work in this chapter was originally published in (149). 104 4.2 Method In this section we consider the NLSE with distributed coefficients in (3+1)D, with a linear potential of the form (4): i∂tu+ β(t) 2 ∆u+ χ(t)|u|2u+ ǫ(t)(x+ y + z)u = iγ(t)u. (4.1) Here t is time, ∆ = ∂2x+∂ 2 y+∂ 2 z is the 3D Laplacian, r = √ x2 + y2 + z2 is the position coordinate, and ǫ(t) stands for the strength of the linear potential as a function of time. The strength of the linear potential is related to the strength of the electric field established in the medium. Hence, we will often use the term “electric field” to refer to ǫ. The direction of the electric field is purposely chosen to symmetrize the three coordinates, though other choices are possible. We assume ǫ(t) 6= 0, otherwise the solutions reduce to those in Sec. 2.3.1. The functions β(t), χ(t), and γ(t) stand for the diffraction, nonlinearity, and gain/loss coefficients, respectively. As in the previous chapters, we search for the solution of Eq. (4.1) in the form (145): u(x, y, z, t) = A(x, y, z, t) exp (iB(x, y, z, t)), (4.2) where: A = f1(t)F (θ) + f2(t)F −1(θ), (4.3) θ = k(t)x+ l(t)y +m(t)z + ω(t), (4.4) B = a(t)(x2 + y2 + z2) + b(t)(x+ y + z) + e(t). (4.5) As in the previous chapters, f1, f2, k, l, m, ω, a, b, e are parameter functions to be determined, and F is one of the Jacobi elliptic functions (JEFs). When Eqs. (4.3)-(4.5) are plugged into Eq. (4.1) and the F-expansion technique (FET) is applied (143), the following set of differential equations for the parameter functions 105 is obtained: df1 dt + 3aβf1 − γf1 = 0, (4.6) df2 dt + 3aβf2 − γf2 = 0, (4.7) dk dt + 2kaβ = 0, (4.8) dl dt + 2laβ = 0, (4.9) dm dt + 2maβ = 0, (4.10) da dt + 2βa2 = 0, (4.11) db dt + 2βab− ǫ = 0, (4.12) dω dt + β(k + l +m)b = 0, (4.13) de dt + β 2 ( 3b2 − (k2 + l2 +m2)c2 )− 3χf1f2 = 0. (4.14) In addition, we obtain two algebraic equations connecting functions f1 and f2 with the coefficients of the Jacobi equation: f1 ( β(k2 + l2 +m2)c4 + χf 2 1 ) = 0, (4.15) f2 ( β(k2 + l2 +m2)c0 + χf 2 2 ) = 0. (4.16) The major change from the system in Sec. 3.3.1 is that equations (4.11) and (4.12) are different. These differences will produce a significant change in the form of solutions. Additionally, one should remark that for b0 = 0 the waves travel in a straight line for both constant β and β = β0 cosΩt. For constant β and b0 = 1, the solution is linear as a function of the longitudinal variable, unlike the case β = β0 cosΩt depicted in Fig. (2.5) where the solution oscillates sinusoidally. 106 4.3 Results Since Eqs. (4.6)-(4.10) and Eqs. (4.15)-(4.16) are the same as in (142), i.e. Sec. 2.3.1, the solutions to all parameters except for b, ω and e are the same: f1(t) = (α) 3/2f0 exp ∫ z 0 γdz, (4.17) f2(t) = δ √ c0 c4 f1, (4.18) a(t) = αa0, (4.19) k(t) = αk0, (4.20) l(t) = αl0, (4.21) m(t) = αm0, (4.22) where α = (1 + 2a0 ∫ t 0 βdt)−1 is the normalized chirp function and δ = 0,±1 and the subscript 0 denotes the value of the given function at t = 0. The integrability condition for χ is also the same: χ = −βc4(k20 + l20 +m20)f−20 e−2 ∫ t 0 γdt/α. (4.23) Equations (4.17)-(4.23) hold for arbitrary β(t) and ǫ(t). We now proceed to solve the remaining equations for four distinct cases: constant ǫ and β; constant ǫ, β = β0 cosΩt; constant β, ǫ = ǫ0 cosΩt; and β = β0 cosΩt, ǫ = ǫ0 cosΩt. We note that in all four subsequent cases the solutions reduce to those found in Sec. 2.3.1 for ǫ = 0. 4.3.1 Case 1: Constant ǫ and β For the case of constant ǫ and β we obtain: b = α(b0 + ǫt+ a0βǫt 2), (4.24) ω = ω0 − αβ(k0 + l0 +m0) ( b0t+ ǫt2 2 ) , (4.25) e = e0 + αβt 2 ( q − ( 3b20 + 3b0ǫt+ ǫ 2t2 + a0βǫ 2t3 2 )) , (4.26) where the normalized chirp function is now α = (1 + 2a0βt) −1 and q = (k20 + l 2 0 + m20)(c2 − 6δ √ c0c4). 107 In Fig. (4.1) (a), (b) we see the effect of adding the electric field on a solitary wave, when the elliptic modulus M of JEFs equals 1. The solution moves away parabolically from the center. The effect of chirp can be seen in Fig. (4.1) (c), (d). As a consequence of the chirp the solution decays. The same will be true for the third case: constant β and ǫ = ǫ0 cosΩt. Both of these effects can be mitigated by setting the value of γ to be γ(t) = 3/(2t). In this way the expression for the chirp will cancel out in the formula for f and g. This idea is first used in (145), and is in detail described in this Thesis in Sec. 3.3. The value to which the strength of the signal converges is 1/(2a0β) (3/2). Dark and bright soliton solutions with chirp and gain are given in Fig. (4.1) (e), (f). Naturally, when gain is present, the power will increase with time. In Fig. (4.2) we see the same solutions as in Fig. (4.1) except that traveling waves are shown, i.e. M = 0.995 < 1. We see by comparing Fig. (4.2) (a), (b) with Fig. (4.2) (c), (d) that the effect of the chirp is not only to reduce intensity but also to broaden the waves. In Fig. (4.2) (e), (f) we see once again the effect of adding gain. Note that the picture is not symmetrical with respect to the transverse variable, but that the central crest curls to the right (towards positive values of the transverse variable), and along with it all the other crests. 4.3.2 Case 2: Constant ǫ and β = β0 cosΩt For the case of constant ǫ and β = β0 cosΩt we obtain: b = α ( b0 + ǫt+ 2a0β0ǫ Ω2 (1− cosΩt) ) , (4.27) ω = ω0 + αβ0(k0 + l0 +m0) · (4.28) (ǫ(1− cosΩt)− Ω(b0 + ǫt) sinΩt) /Ω2, e = e0 + αβ0 2Ω4 ( (((q − 3b20)Ω3 + 3ǫ2Ω + 12a0β0ǫ2Ωt+ 6b0ǫΩ3 (4.29) −3ǫ2Ω2t2) sinΩt+ 6ǫ2Ω2t− 12a0β0ǫ2 − 6b0ǫΩ2) cosΩt −3a0β0ǫ2 cos(2Ωt) + 15a0β0ǫ2 + 6b0ǫΩ2 ) , where α = (1 + 2a0β0 sinΩt/Ω) −1 is the normalized chirp function. In Case 2, we see the effect of a sinusoidal form of β, analyzed in great detail in Sec. 2.3.1, in the presence of a constant electric field. In Fig. (4.3) (a), (b) 108 Figure 4.1: :Soliton solutions for β and ǫ constant as functions of time. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) is presented as a function of k0x+ l0y +m0z and t. Other parameters are: M = 1, β = 1, b0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0.2, δ = 0. Figures (a) and (b) are with no chirp or gain: a0 = 0, γ(t) = 0, Figures (c) and (d) with chirp: a0 = 0.1, γ(t) = 0 and Figures (e) and (f) are with chirp and gain: a0 = 0.5, γ(t) = 3/(2t). 109 Figure 4.2: Traveling wave solutions for β and ǫ constant as functions of time. The parameters are the same as in Fig. (4.1) except for M = 0.995. 110 we see that the effect of the electric field is to widen the amplitude of the soliton. This widening effect happens even in the absence of b0. For b0 = 0 without the electric field, the solution obtained is that of a signal moving in a straight line. In Fig. (4.3) (c), (d) we see the effect of adding chirp, which roughly corresponds to what was obtained in (142). Since β has the form of a sine function the chirp does not cause the signal to decay, given sufficiently low a0 to avoid singularities, but rather the amplitude of the signal oscillates. As in Sec. 2.3.1 there is the effect of stretching along the transverse variable, i.e. the solution is no longer periodic with respect to the transverse variable. In Fig. (4.3) (e), (f) we see the same effect for the corresponding traveling wave. 4.3.3 Case 3: Constant β and ǫ = ǫ0 cosΩt For the case of constant β and ǫ = ǫ0 cosΩt we obtain: b = α ( b0 + ǫ0 sinΩt Ω (1 + 2a0βt) + 2a0βǫ0 Ω2 (cosΩt− 1) ) , (4.30) ω = ω0 − αβ(k0 + l0 +m0) ( b0t+ ǫ(1− cosΩt)/Ω2 ) , (4.31) e = e0 + αβ0 2Ω4 ( (q − 3b20)Ω4t+ 36a0βǫ20 − 24β0ǫ0Ω2 − 6ǫ20Ω2t (4.32) −12a0βǫ20Ω2 + (24b0ǫ0Ω2 − 48a0βǫ20) cosΩt+ 12a0βǫ20 cos(2Ωt) +(3ǫ20Ω + 6a0βǫ 2 0Ωt) sin(2Ωt) ) , where α = (1 + 2a0βt) −1 is the appropriate normalized chirp function. We now analyze the results for Case 3. In the absence of b0 the solutions oscillate as a sine wave, as in Fig. (4.4) (a) and (b). However, if b0 6= 0, then the solution veers toward positive values of the transverse variable, as in Fig. (4.4) (c) and (d). In Fig. (4.4) (e) and (f) we see the traveling wave solutions for b0 6= 0. The presence of chirp causes the solutions to decay, as in the first case. However, the same trick can be applied as in Case 1: adding the right amount of gain to counteract the loss of signal due to chirp. We see this done in Fig. (4.5). 111 Figure 4.3: Soliton and traveling wave solutions for β = β0 cosΩt and ǫ constant. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) presented as a function of k0x+ l0y+m0z and t. Coefficients: β0 = 1, γ(t) = 0, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ = 0.1, Ω = 1, γ(t) = 0, δ = 0. Figures (a), (b) depict solitary waves with no chirp: a0 = 0, M = 1, Figs. (c), (d) solitary waves with chirp: a0 = 0.1, M = 1, and Figs. (e), (f) traveling waves with chirp: a0 = 0.1, M = 0.99999. 112 Figure 4.4: Soliton and traveling wave solutions for β constant and ǫ = ǫ0 cosΩt. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) presented as a function of k0x + l0y +m0z and t. Other parameters are: β = 1, γ(t) = 0, a0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ0 = 1, Ω = 1, δ = 0. Figures (a), (b) depict solitary waves with no b0: M = 1, b0 = 0, Figures (c), (d) solitary waves with b0: M = 1, b0 = 1, and Figures (e), (f) traveling waves with b0: M = 0.99, b0 = 1. 113 Figure 4.5: Soliton and traveling wave solutions for β constant and ǫ = ǫ0 cosΩt with chirp and gain. The parameters are the same as in Fig. (4.4) except for a0 = 0.5, γ(t) = 3/(2t). 114 4.3.4 Case 4: β = β0 cosΩt and ǫ = ǫ0 cosΩt Finally, for the case of β = β0 cosΩt and ǫ = ǫ0 cosΩt we obtain: b = α ( b0 + ǫ0 Ω sinΩt+ a0β0ǫ0 Ω2 (sin2Ωt) ) , (4.33) ω = ω0 − αβ0(k0 + l0 +m0)sinΩt Ω ( b0 + ǫ0 2Ω sinΩt ) , (4.34) e = e0 + αβ0 sinΩt 16Ω4 (8(q − 3b20)Ω3 − 4ǫ20Ω− 3ǫ0(a0β0ǫ0 + 8b0Ω2) sinΩt(4.35) +a0β0ǫ 2 0 sin(3Ωt)), where α = (1 + 2a0β0 sinΩt/Ω) −1 again is the normalized chirp function. In this case the oscillation of β and ǫ combine, as shown in Fig. (4.6). In Fig. (4.6) (c) and (d) we see the effect of adding chirp. As can be seen, chirp warps the amplitude of the signal in surprising ways. Finally, in Fig. (4.6) (e) and (f) we see the same results as in Fig. (4.6) (c) and (d) but for the traveling waves. 115 Figure 4.6: Soliton and traveling wave solutions for β = β0 cosΩt and ǫ = ǫ0 cosΩt as functions of time. Intensity |u|2 for F = sn in Figures (a), (c) and (e) (F = tanh for M = 1) and F = cn in Figures (b), (d) and (f) (F = sech for M = 1) presented as a function of k0x+ l0y+m0z and t. Coefficients: β0 = 1, γ(t) = 0, a0 = 0, b0 = 1, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, ǫ0 = 3, Ω = 1, δ = 0. Figures (a), (b) depict solitary waves with no chirp: a0 = 0, M = 1, Figures (c) and (d) depict solitary waves with chirp: a0 = 0.1, M = 1 and Figures (e) and (f) depict traveling waves with chirp: a0 = 0.1, M = 0.999 116 Chapter 5 Two-component systems 5.1 Introduction In this section the results on two-component systems will be covered (150). Of particular interest is to extend the validity of these exact solutions to a wide range of systems. A system of great interest is the case of two co- and counter-propagating beams coupling with each other through the Kerr nonlinearity (192). The interaction of two beams produces various forms of instability and bifurcation into more chaotic regimes. Also, the case of linearly coupled beams has been extensively studied by Malomed and cases of spontaneous symmetry breaking have been uncovered. In this Thesis, however, we will only concern ourselves with finding novel exact solutions, while the stability of these solutions will be the topic of subsequent research. For simplicity’s sake, we consider only the time-independent (1+1)D case. 5.2 Nonlinear Coupling For the case of two linearly coupled beams, we consider the standard two component nonlinear Schro¨dinger equation with cross-phase modulation in (1+1)-D: i∂zu1 + β(z) 2 ∂2xu1 + χ(z) (|u1|2 + c|u2|2)u1 = iγ(z)u1, (5.1) si∂zu2 + β(z) 2 ∂2xu2 + χ(z) (|u2|2 + c|u1|2)u2 = siγ(z)u2, (5.2) 117 which describes a system of two interacting light beams, u1 and u2, in a medium with Kerr nonlinearity. The functions, as in the previous chapters, β, χ, and γ stand for the diffraction/dispersion, nonlinearity, and gain coefficients, respectively. The coefficient s determines whether the two beams are co-propagating in which case s = 1 or counter-propagating in which case s = −1. It is not to be confused with the coefficient determining whether the dispersion is normal or anomalous used in other chapters. As in the previous chapters, we will separate each of the two beams into the amplitude and phase. 5.2.1 Application of the F-expansion technique We will assume a general ansatz of the following form: u1(z, x) = A1(z, x) exp (iB1(z, x)), (5.3) u2(z, x) = A2(z, x) exp (iB2(z, x)). (5.4) The forms of A1, A2, B1 and B2 are assumed to be: A1 = f1(z)F (θ) + f2(z)G(θ), (5.5) A2 = g1(z)G(θ) + g2(z)F (θ), (5.6) θ = k(z)x+ ω(z), (5.7) B1 = a1(z)x 2 + b1(z)x+ e1(z), (5.8) B2 = a2(z)x 2 + b2(z)x+ e2(z), (5.9) where F andG are two suitable Jacobi elliptic functions, to be determined, satisfying the following differential equations:( dF dθ )2 = c0 + c2F 2 + c4F 4, (5.10)( dG dθ )2 = d0 + d2G 2 + d4G 4. (5.11) By applying the F-expansion method and using the principle of harmonic bal- ance as described in Sec. 2.3.1, we obtain the following equations for f1, f2, g1, g2, 118 ω and k: dfj dz + a1βfj − γfj = 0, (5.12) s dgj dz + a2βgj − sγgj = 0, (5.13) dk dz + 2ka1β = 0, (5.14) s dk dz + 2ka2β = 0, (5.15) dω dz + βkb1 = 0, (5.16) s dω dz + βkb2 = 0, (5.17) where j = 1, 2. The key difference is that we now treat the equations as polynomials in two JEFs, F and G, instead of one. From Equations (5.14) and (5.15) one obtains a1 = sa2 = a and from (5.16) and (5.17) one obtains b1 = sb2 = b. Hence, we ultimately have just two equations for the four parameters a1, a2, b1 and b2: da dz + 2βa2 = 0, (5.18) db dz + 2βab = 0. (5.19) Note that the chirp functions α1 = ( 1 + 2a10 ∫ z 0 βdz )−1 and α2 = ( 1 + 2sa20 ∫ z 0 βdz )−1 are identical for the two parameters a1 and a2, i.e. α1 = α2 = α. We thus obtain the solutions for fj, gj, k, ω, aj and bj, where j = 1, 2: f1(z) = (α) 1/2f10 exp (∫ z 0 γdz ) , (5.20) f1(z) = (α) 1/2f20 exp (∫ z 0 γdz ) , (5.21) g1(z) = (α) 1/2g20 exp (∫ z 0 γdz ) , (5.22) g2(z) = (α) 1/2g20 exp (∫ z 0 γdz ) , (5.23) k(z) = αk0, (5.24) ω(z) = ω0 − αk0b0 ∫ z 0 βdz, (5.25) a1(z) = sa2(z) = a(z) = αa0, (5.26) b1(z) = sb2(z) = b(z) = αb0, (5.27) 119 where j = 1, 2 and α1 = ( 1 + 2a10 ∫ z 0 βdz )−1 is the chirp function. When f2 = g1 = g2 = 0 Eqs. (5.20)-(5.27) reduce to Eqs. (2.19)-(2.24) for the case of a single one-dimensional beam (N = 1) with ǫ = 0. Obviously, the work here is not complete. One needs to find the relationship between f10, f20, g10, and g20, as well as he formula for e. This is determined from equations which are analogous to Equations (2.16)-(2.18). In the case of the co- and counter-propagating beams these equations will induce constraints on the forms of F and G that we can use, as well as the associated parameters c0, c2, c4, d0, d2 and d4. For the equation analogous to Eq. (2.16) we obtain: de1 dz + β 2 ( b2 − kc2 ) = 0, (5.28) de1 dz + β 2 ( b2 − kd2 ) = 0, (5.29) s de2 dz + β 2 ( b2 − kc2 ) = 0, (5.30) s de2 dz + β 2 ( b2 − kd2 ) = 0, (5.31) from which it follows that c2 = d2. Ultimately, one obtains: e1(z) = e10 + 1 2 (c2k 2 0 − b20)α ∫ z 0 βdz, (5.32) e2(z) = e20 + s 1 2 (c2k 2 0 − b20)α ∫ z 0 βdz. (5.33) The equations analogous to Eqs. (2.17)-(2.18) are: f1 ( βc4k 2 + χf 21 + cχg 2 2 ) = 0, (5.34) f2 ( βd4k 2 + χf 22 + cχg 2 1 ) = 0, (5.35) g1 ( βd4k 2 + χg21 + cχf 2 2 ) = 0, (5.36) g2 ( βc4k 2 + χg22 + cχf 2 1 ) = 0. (5.37) Also, an additional set of constraints emerges: 3χf 21 f2 + 2cχf1g1g2 + cχf2g 2 2 = 0, (5.38) 3χf 22 f1 + 2cχf2g1g2 + cχf1g 2 1 = 0, (5.39) 3χg21g2 + 2cχg1f1f2 + cχg2f 2 2 = 0, (5.40) 3χg22g1 + 2cχg2f1f2 + cχg1f 2 1 = 0. (5.41) 120 From Eqs. (5.34)-(5.37) one obtains: (c− 1)f 21 = (c− 1)g22, (5.42) (c− 1)f 22 = (c− 1)g21, (5.43) from which: f1 = ±g2 and f2 = ±g1, or c = 1. We analyze these two cases separately: 5.2.2 Case 1: f1 = ±g2 and f2 = ±g1. Let us assume g2 = δf1 and f2 = ǫg1, where ǫ, δ = ±1. We then obtain from Eqs. (5.38)-(5.41): (3 + c)ǫ+ 2cδ = 0 (5.44) (3 + c)δ + 2cǫ = 0 (5.45) 3 + c+ 2cǫδ = 0. (5.46) The only solution to this system of equations given the initial constraints on ǫ and δ are: ǫ = −δ, (5.47) c = 3. (5.48) Hence, we have ultimately obtained spatio-temporal travelling wave solutions, but only for c = 3. An additional constraint also emerges: f1 g1 = φ √ c4 d4 , (5.49) where φ = ±1. If c4 = d4, i.e. combined with c2 = d2 we have F = G, then one obtains either f1 = −f2, or g1 = −g2, i.e. either A1 = 0 or A2 = 0, which merely reduces to the (1+1)D NLSE. Hence, suitable pairs of JEFs F and G must satisfy c2 = d2, c4 6= d4 and in addition c4d4 > 0 for the ratio in the square root appearing in Eq. (5.49) to be positive. The combinations of allowable forms of F and G (excluding those that are symmetrical) are given in Table (5.1). In cases 1, 2, 3 and 4 we must have M 6= 1, while in cases 5 and 6 we must have M 6= 0. 121 Table 5.1: Possible choices of F and G F G 1 sn ns 2 sn dc 3 cd ns 4 cd dc 5 dn nd 6 sc cs Finally, the integrability condition on χ is found to be: χ = −βαc4k 2 0 exp (−2 ∫ z 0 γdz ) (c+ 1)f 210 . (5.50) 5.2.3 Case 2: c = 1. Let us now, on the other hand, assume c = 1. Setting g2 = δf1 and f2 = ǫg1, where ǫ and δ are real parameters not necessarily equal to ±1, one obtains the following equations: 3ǫ+ 2δ + ǫδ2 = 0, (5.51) 3δ + 2ǫ+ ǫ2δ = 0, (5.52) 3ǫ2 + 2ǫδ + 1 = 0, (5.53) 3δ2 + 2ǫδ + 1 = 0. (5.54) By equating the right-hand side (RHS) of Eq. (5.53) with the RHS of Eq. (5.54), one obtains δ = ±ǫ, and then from Eqs. (5.53)-(5.54) one obtains either 5ǫ2 + 1 or ǫ2 + 1 = 0. Since neither of these two equations has a real solution, it follows that solutions of the form given in Eqs. (5.5)-(5.9) are not possible for c = 1. Hence, the present method does not provide solutions for the Manakov system of equations. 5.2.4 Results We have thus found from our analysis in Sec. 5.2.1 that only for c = 3 one obtains the traveling wave solutions of the form given in Equations (5.5)-(5.9) by the present 122 method. This, of course, does not preclude the existence of solutions for other values of c. For example, the value of c = 1 corresponds to the simple Manakov model. It is known that the Manakov model is integrable by the inverse scattering method (54). As for the feasibility of constructing materials for which c = 3, in (52) it has been obtained that by using periodically poled photorefractive media one can eliminate SPM, and since this process can be graded a full range of real values of c can be achieved. We now present the traveling wave solutions obtained. Of most interest are solutions in case 5 since for M 6= 1 functions dn and nd do not have singularities. We present the results in Fig. (5.1). The general form of the solutions is the same for both co- and counter-propagating beams. We see that the solution remain periodic and do not decay in the absence of loss. Figure 5.1: Traveling wave solutions for F = dn and G = nd constant as functions of time. Intensities: (a) |u1|2 and (b) |u2|2 are presented as a function of k0x and z for β(z) = β0 cosΩz. Other parameters are: M = 0.99, b0 = 0, e0 = 0, k0 = l0 = m0 = 1, ω0 = 0, Ω = 1, β0 = 1, γ = 0, ǫ = φ = 1 and δ = −1. In Fig. (5.2) we see the effects of chirp on the solutions. The effects are similar to those described in (142) for one-component system, namely a modulation of amplitude and deformation along the k0x axis. Finally, in Fig. (5.3) we see the combination of two different types of functions giving us novel solutions to the coupled NLSE, albeit with a singularity. 123 Figure 5.2: Traveling wave solutions as functions of time. The parameters are the same as in Fig. (5.1) except for a0 = 0.1. Figure 5.3: Traveling wave solutions as functions of time. The parameters are the same as in Fig. (5.1) except for F = sn and G = dc. 124 Chapter 6 Stability Analysis 6.1 Introduction In this section the stability of the solutions obtained in the previous chapters will be investigated, most notably the solutions to the Nonlinear Schro¨dinger equation (NLSE) in Sec. 2.3.1 and the solutions to the Gross-Pitaevskii (GP) equations in Sec. 3.3. The stablity of exact soliton solutions to the NLSE is a very important question to be answered (193). In (1+1) dimensions, bright and dark soliton solutions of the NLSE in Kerr medium with cubic nonlinearity are unconditionally stable for the, respectively, self-focusing and self-defocusing nonlinearity (4). However, in homoge- nous bulk media with a self-focusing cubic Kerr nonlinearity one cannot have un- conditionally stable two and three dimensional solutions of the NLSE. Nevertheless, great interest has been generated when it was suggested that the two dimensional generalized NLSE with varying coefficients may lead to stable 2D solitons (80). The stabilizing mechanism has been the sign-alternating Kerr nonlinearity in a layered medium. In this Thesis, the solutions obtained frequently utilized relationships for the nonlinearity term χ that resulted in the alternating sign of the Kerr nonlinearity (2.27) and therefore it is worth investigating whether the solutions obtained in this thesis are stable. Initial work on stability analysis was done in (142) and the results are shown 125 Figure 6.1: Numerical simulation of the light bullet from Fig. (2.5)(b). Initial data from Eq. (2.56) are propagated according to Eq. (2.30) for 90 diffraction lengths along the z axis. Only the dependence on t is shown. The initial profile is noted by open circles. The curves to the left present intensity profiles at the left turning point, the curves to the right the profiles at the right turning point. The curves at the center are snapshots of the profiles passing approximately through the point t = 0 (i.e., the frames closest to t = 0 are recorded). Three sets of 15 profiles are overlapped at different z points, to show that no instabilities develop. in Fig. (6.1). We showed that the bright soliton solutions remain stable for a very long time. The solutions for the bright soliton were run in a computer simulation of the propagation according to Eq. (2.30) and showed no sign of blowing up after 90 diffraction lengths. However, numerical solutions are limited in their power to conclusively determine the stability of a solution, therefore analytical methods must be used to ascertain the stability of the obtained solutions. We will use standard techniques to explore the modulational instability of the solitary wave solutions obtained in Sec. 2.3.1 and Sec. 3.3, using a variational approach. Localized two- and three-dimensional solutions of the cubic nonlinear Schro¨dinger equation in the works (143; 142; 144) are extended one dimensional solutions. Namely, the intensity |u| of solutions from the are homogenous in one (or two) out of two (or three) dimensions. It is in this direction (or plane), due to the nonlinearity, that modulation instability can develop. For this reason, particular interest is the analysis of modulation (in)stability of solution in the direction of 126 homogeneity of |u|. 6.1.1 Modulational Instability In this subsection we will develop the basic theory of modulational instability (MI). The main idea of the study of MI is to introduce perturbations in the solutions of the given partial differential equation (PDE), plug in the perturbed solution into the PDE, keeping only the first order terms, and determine the spatial and temporal evolution of the said perturbation. If the amplitude of the perturbation grows exponentially, or even linearly, then this is an indication of MI. Conversely, in the case of modulational stability (MS), the amplitude of the perturbation will typically oscillate periodically. Even if the amplitude reaches a very high value, as long as it is periodical, the solution to the PDE is said have MS. While other types of instability are not necessarily excluded, modulational instability remains a powerful tool for analyzing PDEs, as well as ascertaining the feasibility of finding the given class of solutions. As an example, let us consider the nonlinear Schro¨dinger equation (NLSE) of the following form: i ∂E0 ∂z + β(z) 2 ( ∂2E0 ∂x2 + ∂2E0 ∂y2 + s ∂2E0 ∂t2 ) + χ(z) |E0|2E0 = 0, (6.1) where E0 is the unperturbed solution to Eq. (6.1). Adding now a perturbation to the solution v, we now plug in the perturbed solution: E = E0 + v, (6.2) into the equation and keeping only the first order terms in v, we obtain the following equation: i ∂v ∂z + β(z) 2 ( ∂2v ∂x2 + ∂2v ∂y2 + s ∂2v ∂t2 ) + χ(z) ( 2|E0|2v + E20v∗ ) = 0, (6.3) where v∗ is the complex conjugate of v. Equation (6.3) is typically split into real and imaginary parts by assuming v = vr + ivi. The obtained system of equations becomes: −∂vi ∂z + β(z) 2 ( ∂2vr ∂x2 + ∂2vr ∂y2 + s ∂2vr ∂t2 ) + χ(z) ( 2|E0|2 + E20 ) vr = 0, (6.4) ∂vr ∂z + β(z) 2 ( ∂2vi ∂x2 + ∂2vi ∂y2 + s ∂2vi ∂t2 ) + χ(z) ( 2|E0|2 − E20 ) vi = 0. (6.5) 127 If exact solutions cannot be found to the coupled system of differential equations, numerical techniques are typically employed. 6.2 Nonlinear Schro¨dinger equation and its trans- formation We confine our analysis to the (2+1) and (3+1)-dimensional NLSEs considered in (143; 142; 144), and use the notation introduced there and in this Thesis. We consider the generalized NLSE with Kerr nonlinearity developed in Secs. 2.3.1-2.3.4: i ∂u ∂z + β(z) 2 ( ∂2u ∂x2 + ∂2u ∂y2 + s ∂2u ∂t2 ) + χ(z) |u|2 u = iγ(z)u. (6.6) Our goal will be to verify whether the solutions to Eq. (6.6) developed in Secs. 2.3.1-2.3.4 are MS or with MI. Here, in Eq. (6.6), s = −1 for the normal dispersion and s = 1 for the anomalous dispersion. For s = 0 we have the two-dimensional time independent NSLE which will also be covered in our stability analysis. The coefficient of nonlinearity χ is given as χ(z) = −c4β(z) αf 20 χ0 exp(−2 ∫ z 0 γdz), (6.7) where we define: χ0 = (k 2 0 + l 2 0 + sm 2 0) (6.8) and α is the chirp function as defined in Eq. (2.26). The parameters a, b, k, l, and m are as defined in Eqs. (2.45)-(2.53) and Eqs. (2.68)-(2.76) for the anomalous and normal dispersions, respectively. For convenience, we will rewrite the parameter e as follows: e = (α/2) ( c− (2 + s)b20 ) ∫ z 0 βdz = e0 + αc 2 ∫ z 0 βdz = e0 + c 2 ∫ z 0 βα2dz, (6.9) with e0 defined accordingly (not as the starting value which we have assumed to be 0), and where we define c = c2 − 6ǫ√c0c4. We will limit our attention to the cases where ǫ = 0, hence c = c2. We now make the following substitution: G = u exp ( − ∫ z 0 γdz ) exp (−i(a(x2 + y2 + st2) + b(x+ y + t) + e0))/(f0α3/2) (6.10) 128 to reduce the equation Eq. (6.6) to a more manageable form. In addition, we perform a gauge transformation on the coordinates: x → x′ = α(x− ς), (6.11) y → y′ = α(y − ς), (6.12) t → t′ = α(t− sς), (6.13) z → z′ = ∫ z 0 αβ2dz, (6.14) where ς(z) = b0 ∫ z 0 βdz. By plugging in all the new variables into Eq. (6.6) we obtain the following equation: i ∂G ∂z′ + 1 2 ( ∂2G ∂x′2 + ∂2G ∂y′2 + s ∂2G ∂t′2 ) − c4χ0 |G|2G = 0. (6.15) Equation (6.15) is much more suitable for our stability analysis than Eq. (6.6) because all the coefficients next to G and its derivatives are constant. Also, the wave propagation now necessarily happens on a straight line, unlike in solutions shown in (142; 144). The stationary solutions: F = G exp ( −ic2 ∫ z 0 βα2dz/2 ) , (6.16) where G is the solution of Eq. (6.15), contain the whole range of solutions from (142) and (144), for different values of c2 and c4, in the direction (k0, l0,m0). Since ǫ = 0, F will be equal to some JEF. Without loss of generality, we can put f0 = 1 and also k0 = 0, l0 = 0, m0 = 1, for temporal and k0 = 0, l0 = 1, m0 = 0, for spatial solutions. Adjusting parameters (k0, l0,m0) to these values corresponds to the rotation and re-scaling of the coordinate system (x′, y′, t′). We note that in this case, χ0 = 1 for spatial solitons, whereas χ0 = s for temporal solitons. For both the spatial and temporal solitons x′ will be chosen to be the direction of modulation, as a perpendicular direction is needed for study (any will suffice) and x′ is perpendicular to both, given the conventions we have chosen. Also, if the soliton is spatial, perturbations along the temporal axis t′ are also possible. In that case, to unify our results, we will consider the x′ and t′ to be swapped, so that x′ always denotes our axis of perturbation. 129 6.3 Variational approach of the modulation sta- bility We will now use the variational approach to examine MI, following (194). We will analyze two periodic planar solutions: F = sn (dark soliton) F = cn (bright soliton) which reduce in the limit of M = 1 to tanh and cosh−1, respectively. We will study spatial and temporal periodic solutions for both cases: the normal and the anomalous dispersion. The idea of the variational approach is to introduce a perturbation in our solu- tions. In other words, we assume: G = G0 (1 + (Ur(z) + iUi(z)) cos (Kx)) , (6.17) where G0 is the unperturbed solution of Eq. (6.15), K is the wave number of the perturbation and x is the direction in which the perturbation occurs, as described in the previous section. We then construct, according to standard procedure, the Lagrangian corresponding to Equation (6.15): L = i 2 ( G ∂G∗ ∂z′ −G∗∂G ∂z′ ) + 1 2 |∇G|2 + c4χ0 |G|4 , (6.18) where G∗ is the complex conjugate of G, and |∇G|2 = |∂G/∂x′|2 + |∂G/∂y′|2 + s |∂G/∂t′|2. We now perform an averaging of the Lagrangian over all three transverse coordinates to obtain: 〈L〉 = ∫ ( i 2 ( G ∂G∗ ∂z′ −G∗∂G ∂z′ ) + 1 2 |∇G|2 + c4χ0 |G|4 ) dx′ dy′ dt′. (6.19) Note that the Lagrangian has been averaged over one period of perturbation in the direction of perturbation and for all other transverse directions it has been averaged from −2K(M) to 2K(M) for F = cn and from −K(M) to 3K(M) for F = sn (these boundaries converge to −∞ and ∞ for M = 1, i.e. for the solitary waves). Here K(M) = ∫ pi/2 0 (1 −M sin2 t)−1/2dt is the complete elliptic integral of the first kind. The total action is now defined as: Λ = ∫ ∞ −∞ 〈 L〉dz. (6.20) It remains invariant in the transformation of coordinates from u to G. 130 Substituting the new formula for G into Eq. (6.60) into the effective Lagrangian given in Eq. (6.20), we vary Ur and Ui in the standard procedure (194) to obtain the Euler-Lagrange ordinary differential equations for Ur = Ur(z) and Ui = Ui(z), as follows: ∂ β∂z (αUr) = 1 2 K2 (ασUi) , (6.21) ∂ β∂z (ασUi) = −1 2 ( K2 − σσηdα2 ) (αUr) , (6.22) where η = −χoc4, ση = sgn(η), d = d(M)cn = 8 3 (2M − 1)(E(M)− E(am(5K(M)|M)|M))−2(2− 5M + 3M2)K(M) E(M)− E(am(5K(M)|M)|M))− 4(M − 1)K(M)) , (6.23) for F = cn(·|M), and d = d(M)sn = 8 3 (M + 1)E(am(4K(M)|M)|M)− 2(2 +M)K(M) (E(am(4K(M)|M)|M)− 4K(M)) , (6.24) for F = sn(·|M). Here, M is the parameter of the JEF, K(M) = F (π/2|M) and E(M) = E(π/2|M) are the complete elliptic integrals of the first and second kind, respectively, F (u|M) = ∫ u 0 (1−M sin2 t)−1/2dt and E(u|M) = ∫ u 0 (1−M sin2 t)1/2dt are the incomplete elliptic integrals of the first and second kind, respectively, and am(u,M) = F−1(u,M) is the amplitude of Jacobi elliptic functions. The parameter σ is defined as follows: σ = 1 for the perturbations along the spatial coordinates (spatial perturbations) and σ = s for the temporal perturbation of spatial solitons. The dependence of coefficients d (m) cn and d (m) sn , on the elliptic modulus of the JEF M (0 < M < 1) is shown in Fig. (6.2)(a). For bright solitons d (M=1) cn = 8/3 and for dark solitons d (M=1) sn = 4. Equations (6.21)-(6.22) are ODEs which can be solved using standard methods. The solution to Eqs. (6.21)-(6.22) can now be written as: Ur = α −1 (z) ( C1M0,ν(iK 2/2a0α(z)) + C2W0,ν(iK 2/2a0α(z)) ) , (6.25) Ui = 2σ α(z)K2 ∂ β∂z ( α(z)Ur ) , (6.26) where α(z) ≡ (1 + 2a0 ∫ z 0 βdz)−1 ≡ (1 + 2a0ξ)−1 > 0 is the chirp function, M0,ν(z) and W0,ν(z) are Whittaker functions: ν = −1 2 √ 1 + σσηK2d 4a20 , (6.27) 131 Figure 6.2: (a) Nonlinearity parameter d for solutions cn and sn. (b) The growth rate parameter γ for dark an bright solitons as a function of K for the case σση = 1. Modulational instability occurs for values of K depicted on the respective graphs. The solid lines represent the theoretical calculation of K using Eq. (6.33), and the square and circle dots are values of γ measured using numerical simulations, in which the dark and bright solitons, respectively, were perturbed by a small wave of the given wave number K. and constants C1 and C2, in the case of initial conditions Ur(z = 0) = U0 and Ui(z = 0) = 0, become: C1 = −U0 √ π 21+2νΓ(1 + ν) ( W0,ν(iK 2/2a0)− 2 2a0 iK2 W1,ν(iK 2/2a0) ) , (6.28) C2 = U0 √ π 23+2νΓ(2 + ν) ( M0,ν(iK 2/2a0) + 2a0 iK2 (1 + 2ν)M1,ν(iK 2/2a0) ) .(6.29) For small values of ξ = ∫ z 0 βdz << 1, the modulus of the perturbation amplitude can be approximated to within second order of ξ to be: |U | = U0 ( 1 + 2a0ξ − (σσηd/8) ( K2 − σσηd ) ξ2 ) . (6.30) 6.3.1 Case without chirp In the case of without chirp, i.e. a0 = 0, solutions to Eqs. (6.21)-(6.22) become Ur = U0 cosh (γξ) , (6.31) Ui = U0 2γ K2 sinh (γξ) . (6.32) 132 Here, γ is the MI gain: γ = K √ σσηd−K2/2, (6.33) also known as the growth rate parameter, and is not to be confused with with the gain/loss parameter of the NLSE, which will not be used subsequently in this Chapter. The dynamics of the overall evolution of the total perturbation: U = Ur + iUi (6.34) are determined by the the MI gain (4) and the function β(z). For σση = −1 the growth rate parameter γ = iγ¯ = iK √ K2 + d/2 is imaginary for all values of K, and, consequently, the solutions G are modultionally unconditionally stable for any function β(z). This case occurs for temporal cn-solitons with normal dispersion (s = −1) and temporal sn-solitons for s = 1 in the self-defocusing media. For example, if we take the function of β to be β(z) = β0 + β1 sin(2πz/Z), where Z is a constant representing the wavelength of β, then the amplitude of the spatial perturbation satisfies the condition: ( 1 + d/K2 )1/2 ≥ |U/U0| ≥ 1, (6.35) for any values of the parameters β1 and β0 6= 0, and: ( 1 + ( d/K2 ) sin2(γ¯β1Z/π) )1/2 ≥ |U/U0| ≥ 1, (6.36) for β0 = 0. In this case, σση = −1, the solutions are unconditionally stable. In the opposite case, σση = 1, which holds for temporal or spatial cn-solutions for anomalous dispersion s = 1 or temporal sn-solutions if s = −1 in the focusing media, more interesting dynamics of perturbations occur. For the spatial perturbation with wavenumber K > √ d the growth rate is zero (since γ = iγ¯ is imaginary) and the perturbation amplitude has an oscillatory solution in the following range: ( 1− (d/K2) sin2(γ¯β1Z/π))1/2 ≤ |U/U0| ≤ 1. (6.37) If K < √ d, then |U | grows exponentially at a rate γβ0, exhibiting either oscillations (if β1 6= 0) or monotony (if β1 = 0), as seen in Fig. (6.3)(a). Consequently, if β0 6= 0, the solution is unstable. In Fig. (6.2)(b) we plot γ as a function of the wavenumber K for bright (d = 8/3) and dark (d = 4) solitons in the case without 133 Figure 6.3: (a) Theoretical values of log |U/U0| based on Eqs. (6.25)-(6.26) for a0 = 0. (b) Theoretical values of |U/U0| based on Eqs. (6.25)-(6.26) for a0 = 0.1. dispersion management: β(z) = β0 = 1. Numerical simulations of Eq. (6.6) confirm the analitical prediction of the growth rate of γ. Thus, the amplitude of |U | can be made to be stable if the mean value of the management function is zero β0 = 0, and the period of oscillations of β is small (i.e. Z << 1). The variation of the perturbation amplitude in this case is: |U/U0| < ( 1 + ( d/K2 ) sinh2(γ¯β1Z/π) )1/2 . (6.38) The entire stability analysis is presented in Table (6.1). We see that depending on the choice of s, ση and whether the soliton is spatial or temporal, we have eight distinct cases for examining the stability of our solutions. In the case of spatial solitons, perturbations can occur in both the spatial and temporal directions, whereas in the case of temporal solitons they can only occur in the spatial directions. For a soliton to be stable we must have σση = −1 in all directions of perturbation, otherwise it is conditionally stable, i.e. only for β0 = 0. A soliton is dark if the direction of the soliton and the nonlinear term, i.e. ση, are of the opposite sign, otherwise the soliton is bright. In the 2D time-independent case we no longer need to consider temporal solitons and distortions. 134 Table 6.1: Stability cases s ση soliton pert. σ type stability (3D) stability (2D) 1 1 1 S-spatial S T 1 1 cn CS-Conditionally Stable CS 2 1 1 T-temporal S 1 cn CS – 3 1 −1 S S T 1 1 sn S-stable S 4 1 −1 T S 1 sn S – 5 −1 1 S S T 1 −1 cn CS CS 6 −1 1 T S 1 sn CS – 7 −1 −1 S S T 1 −1 sn CS S 8 −1 −1 T S 1 cn S – 6.3.2 Case with chirp In the case where a0 6= 0 we analyze the propagation of Eqs. (6.25)-(6.26) as a function of z, with C1 and C2 defined in Eqs. (6.28)-(6.29). The behavior is similar to the case without chirp and Table (6.1) is also applicable in this case; however, there are a handful of distinctive features. Typical behavior of Eq. (6.34) is presented in Fig. (6.3b), for different values of parameters. If β0 6= 0 we have averaged, i.e. linear, growth of the perturbation amplitude for all values of wave number K. In the case β0 = 0, the perturbation amplitude has oscillatory behavior with a maximum variation that depends on the period Z of management function β. Therefore, if the mean of management function is different from zero, the mod- ulation instability grows linearly and the soliton is unstable. Stabilization can be achieved by reducing the period of the management functions and by satisfying the condition β0 = 0. Finally, we summarize all our findings in Fig. (6.4) for the conditionally stable case. The plots depict the maximum value of the amplitude U in units of U0. The 135 purple region is where U never ascends above its initial value, U0. In the case where In Fig. (6.4)(a) we see that the modes above √ 8/3 ≈ 1.63 are all stable and do not increase the amplitude of our function, while for lower values of K the intensity of the perturbation depends on Ω = 2π/Z, the angular frequency of the oscillation of β. For large values of Ω the intensity of the oscillations is low, while for small values of Ω it blows up, as predicted in our analysis in Sec. 6.3.1. With the addition of chirp we see that the general trends as the magnitude of chirp a0 increases are that the amplitude of oscillations increases for values of K above √ 8/3 and also for higher values of Ω when K < √ 8/3. 6.4 Numerical simulations We will use the split-step fast Fourier transform (SSFFT) method in numerically simulating our solutions. The goal of our computer simulations is to verify the rate of modulational instability growth predicted by the value of parameter γ given in Eq. (6.33). 6.4.1 Split-step fast Fourier method The split-step fast Fourier method was designed to simulate the evolution of PDEs that combine a linear and nonlinear operator (99). It is extensively used in the computer simulations of the NLSE (195; 196) and GP equations (197). The main goal of the split step method is to reduce the accumulation of error in the standard application of the Fast Fourier Trasform (FFT). If we wish to simulate the evolution of a nonlinear partial differential equation along the direction of propagation, we will usually be dealing with both a linear and a nonlinear term in the equation: ∂u ∂z = i(L+N )u, (6.39) where L is a linear operator on u and N is a nonlinear operator on u. Typically, L(u) = β 2 △ for the NLSE equation and L(u) = β 2 △ + V (r) for the GP equation, where △ is the transverse Laplacian and V is a potential function of the transverse 136 Figure 6.4: Maximum amplitude of oscillations plotted against Ω andK. Parameters used: β0 = 0, β1 = 1, d = 8/3, (a) a0 = 0, (b) a0 = 0.02, (c) a0 = 0.1 and (d) a0 = 0.3. The contour boundaries, starting from the color purple are: 1, 2, 5, 10, 20, 50. 137 radius r. In both cases we will have N = χ|u|2. Using the standard operator formalism we can write a formal solution of Eq. (6.39) as: u(z) = e i ∫ z z0 (L+N )dz u(z0). (6.40) Since we will only be taking a small step from z to z+ △ z we can assume the values of the operators to be roughly constant and utilize the Baker-Hausdorf theorem to write: u(z+ △ z) = eiA+iBu(z) ≈ eiB/2eiAeiB/2u(z), (6.41) where A = L △ z and B = N △ z. This approximation is valid up to the third degree of △ z. Higher order approximations can also be found in the literature (198). We now consider the following two differential equations: ∂u ∂z = i 2 Nu = i 2 χ|u|2u, (6.42) ∂u ∂z = iLu. (6.43) We note that the two equations correspond to the application of eiB/2 and eiA respectively for small steps. From Eq. (6.42) and its conjugate it follows that: 0 = ∂u ∂z u∗ + ∂u∗ ∂z u = ∂(|u|2) ∂z . (6.44) In other words, |u| = const., and hence N = const., which significantly simplifies the treatment of Eq. (6.42). As for Eq. (6.43) we apply the FT to transform the equation into the frequency domain: ∂u ∂z = −iβ 2 κ2u, (6.45) where κ = √ κ2x + κ 2 y + κ 2 t is the overall wave number. Thus, combining all these results the new formula for u is equal to: u(z+ △ z) = eiχ|u| 2 △z/2FT−1(e−i β 2 κ2△zFT(eiχ|u| 2 △z/2u(z))), (6.46) where FT(u)(κx, κy, κt, z) = ∫ ∫ ∫ ∞ −∞ u(x, y, t, z)e−(κxx+κyy+κtt)dx dy dt, (6.47) FT−1(u)(x, y, t, z) = 1 (2π)3 ∫ ∫ ∫ ∞ −∞ u(κx, κy, κt, z) · (6.48) e(κxx+κyy+κtt)dκx dκy dκt. 138 In order to apply Eq. (6.46) in a computer simulation we will, following (198), perform a discretization of both the physical space and the momentum space. We will, for convenience restrict ourselves to the interval of 0 to 2π for each spatial and temporal coordinate and separate the physical space in the interval into N discrete pieces, where N is an even number, forming the grid points xi, zi, ti = 2πi/N , where i = 0, 1, . . . , N−1. Hence, Upqr(z) = u(xp, yq, tr, z) are going to be the approximate values of u on the grid points. For the momentum space we will separate the grid into the region from −N/2 to N/2+1, i.e. ki, li,mi = i, where i = −N/2, −N/2+1, . . . , N/2+1, while the Fourier transform will be evaluated as: Upqr(z) = u(kp, lq,mr, z). Taking into account the normalization of the coefficients, the formulas for the Fourier transform now become: Upqr(z) = 1 N3 N−1∑ s,t,v=0 Ustv(z)e −i(kpxs+lqyt+mrtv), (6.49) Upqr(z) = N/2∑ s,t,v=−N/2 U stv(z)e i(ksxp+ltyq+mvtr). (6.50) The computation is significantly speeded up when applying FFT, which can be implemented when N = 2k, for some integer k. The most common algorithm used for computing the FFT is the Cooley-Turkey algorithm (199), originally invented by Gauss, which establishes recurrent relationships between FFTs of composite order N = N1N2, assuming the FFT computation for N1 and N2 is known. The basic formulas developed in this section clearly describe how FFT is imple- mented in simulations of many different nonlinear partial differential equations. 6.4.2 Results We now use Eq. (6.46) to perform computer simulations of the behavior of our solutions when a small perturbation is introduced. Observing the rate of change of the amplitude of G in our simulations, we can then measure the value of γ and compare it with the theoretical expectation given in Eq. (6.33). Our FFT simulations gave us the points on the plot of Fig. (6.2b) which agree well with theoretical expectations denoted by the continuous line. Most importantly, the solutions cease to exponentially increase at precisely the values predicted by the theory of MI. 139 Figure 6.5: Development of modulational instability for the bright soliton for three different values of z. Here, x is the direction of perturbation, y is the direction of the soliton and t is the remaining transverse direction. Lower-wavelength colors (i.e. towards the color red) indicate a higher value of |u|2. Figure 6.6: Development of modulational instability for the dark soliton for three different values of z. Here, x is the direction of perturbation, y is the direction of the soliton and t is the remaining transverse direction. Lower-wavelength colors (i.e. towards the color red) indicate a higher value of |u|2. 140 We see in Figs. (6.5)-(6.6) the main results of our simulations in scenarios involving instability. After initially being unnoticeable (Fig. (a)), we can see the perturbation rapidly increase and ultimately completely abandon the original form of the solution (Figs. (b) and (c)). Figure (6.5) depicts the time evolution of a bright soliton, while Figure (6.6) depicts the time evolution of a dark soliton. Due to the difficulty that arises in the boundary conditions for the dark soliton due to a change of sign of u we have instead run a simulation of two dark solitons with periodic boundary conditions, as is standard practice. 6.5 Analysis of stability for the Gross-Pitaevskii equation In this section we will apply the results and methods of Sec. 6.2 to the solutions to provide a stability analysis of our solutions to the Gross-Pitaevskii (GP) equation obtained in Sec. 3.3. We now examine Eq. (3.1), i.e.: i∂tu+ β(t) 2 ( ∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 ) +χ(t)|u|2u+α(t)(x2+ y2+ z2)u = iγ(t)u. (6.51) The parameter α in this case refers to the strength of the quadratic potential and not the chirp function, as in the previous sections. The main part of our analysis will be to transform the starting Eq. (6.51) into a form more amenable to stability analysis (200). The key differences between Eq. (6.51) and Eq. (6.1), apart from the addition of the quadratic potential, is the change in the longitudinal direction from z to t. Hence, there is no longer a distinction between normal and anomalous dispersion. This greatly simplifies the stability analysis. We will restrict out attention to solitons found in Sec. 3.3.2, as the solutions found in Sec. 3.3.1 do not have a stable amplitude when they are not artificially maintained with a nonzero gain. Formula (6.7) now becomes: χ(z) = −c4β(z) pf 20 χ0 exp(−2 ∫ z 0 γdz), (6.52) where χ20 = k 2 0 + l 2 0 + m 2 0 and p is defined in Eq. (3.29). For convenience we set χ0 = 1, f0 = 1 and β0 = 1, where β0 is the amplitude of β. Our equation for e now 141 becomes: e = (q/2) ( c− 3b20 ) = e0 + qc 2 , (6.53) where q is given in Eq. (3.30) and c = c2−6ǫ√c0c4. Again, for ǫ = 0 we have c = c2. We now define: G = u exp ( − ∫ z 0 γdz ) exp (−i(a(x2 + y2 + st2) + b(x+ y + t) + e0))/(f0p3/2) (6.54) and employ similar changes of coordinates as in Sec. 6.2: x → x′ = p(x− ς) (6.55) y → y′ = p(y − ς) (6.56) t → t′ = p(t− sς) (6.57) z → z′ = ∫ z 0 pβ2dz, (6.58) where p is defined in Eq. (3.29) and ς ′(t) = β(t)(2a(t)ς(t) + b0p(t)). The form of the function ς(t) in the transformation is of no immediate interest, other than the fact that it depends on p(t) and τ(t) which are given in Eq. (3.29) and Eq. (3.32), respectively. The transformation gives us the following equation for G: i ∂G ∂t′ + ( ∂2G ∂x′2 + ∂2G ∂y′2 + ∂2G ∂z′2 ) − c4 |G|2G = 0. (6.59) Qualitatively, this is the same equation as Eq. (6.15), except for the exchange of the longitudinal variable from z to t. The new variable t′, which only depends on t, involves an integral over β that can change sign. This will be important in the analysis of Eq. (6.59). Equation (6.59) is the usual (3+1)D nonlinear Schro¨dinger equation with con- stant coefficients, which is prone to instabilities and the wavefunction collapse. In- stabilities in G translate into instabilities of the general solution u. This would bode disaster for the stability of exact traveling wave and solitary solutions found, were it not for the possibility of diffraction and nonlinearity management (153) in Eq. (6.59), thanks to the form of the primed variables. We find that, for the choice of coefficients α(t) and β(t) made in Sec. 3.3.2, the typical extended soliton solu- tions of Eq. (6.59) do not collapse when perturbed, but keep oscillating in a typical breathing behavior. 142 Without loss of generality we place the z′ axis in the direction of inhomogeneity of our extended solitary solutions. In this manner the z′ variable takes the role of the θ variable. The intensity |G| is homogenous in two of the three spatial dimensions (i.e. in the plane perpendicular to the direction (k0, l0,m0) of inhomogeneity.) It is in this plane that, owing to nonlinearity, the modulational instability can develop. For this reason, of particular interest is the analysis of MI of perturbations in the plane of homogeneity of |G|. We consider the perturbation of G in this plane for the two fundamental solu- tions, the dark F = sn and the bright F = cn solitons, where F = G exp (−iqc2/2) in the form: G = G0 (1 + (Ur + iUi) cos(Kx)) , (6.60) where U(t) = Ur(t) + iUi(t) is the complex amplitude, and K is the wavenumber of the perturbation in the direction perpendicular to z′. In a standard linear stability analysis, as was already done for the NLSE in Sec. 6.2, the perturbation is substi- tuted into Eq. (6.59) and linear first-order differential equations for Ur and Ui are obtained, which are then solved to analytically to yield: Ur(t) = (C1P µ ν (tanh(τ)) + C2Q µ ν (tanh(τ))) p −1, (6.61) Ui(t) = 2 √ 2α0 pK2 ∂ [pUr] ∂τ , (6.62) where P µν and Q µ ν are the associated Legendre functions, with: ν = −1/2 ( 1− √ 1− dK2/2 (α0 − 2a20) ) , (6.63) and µ = iK2/2 √ 2α0. The solutions of Eqs. (6.61) and (6.62) determine the dy- namics of the modulational instability. The constants C1 and C2 are determined by the initial conditions for Ur and Ui, Ur(0) = U0 and Ui(0) = 0. Also, d = −4 for the dark soliton and d = 8/3 for the bright soliton for M = 1. In fact, a similar analysis as in Sec. 6.3.1 can be performed in this case, and one obtains a much more simplified analysis with respect to Table (6.1) since there is no case s = 1, and there are no temporal perturbations. One obtains that the dark solitons are always stable, and the bright solitons are always conditionally stable. We now restrict our attention to the bright solitons. Figure (6.7) depicts a typical evolution of the modulus of the perturbation amplitude |U | for different 143 Figure 6.7: Evolution of the modulus of the perturbation amplitude |U | of bright solitons in time for different values of Ω: a) Ω = 1, b) Ω = 8. In both figures a0 = 0.3. Other parameters: K = 1, α0 = 0.3 [black, (upper) solid line]; K = 1, α0 = 0.1 [black, (upper) dashed line]; K = 4, α0 = 0.3 [red, (lower) solid line] and K = 4, α0 = 0.1 [red, (lower) dashed line]. values of the parameters. In all the cases, we have a periodic dependence in time, with the period 2π/Ω, where Ω is the frequency of the modulation of α and β. For large K we see a superposition of two oscillations, one with the frequency Ω and the other with the frequency ∼ K√K2 − d/2. For small Ω, independent of the value of K, the amplitude of the perturbation may, for a period equal to π/2Ω, grow for several orders of magnitude compared to the initial value. This is expected since the generalized GPE should display sensitivity to low-frequency (long period) perturbations. Such perturbations of the coefficients bring GPE closer to the limit of NLSE with constant coefficients, which naturally is prone to instability and collapse. In contrast to this case, for large Ω the variation of the perturbation amplitude is much smaller. In addition, a decrease in the initial chirp a0 and an increase in the strength of the potential α0 cause a reduction in the variation of the perturbation amplitude. These two trends are clearly seen in Fig. (6.8). In fact, for large Ω ≫ √2α0 and small chirp α0 ≫ 2a20 the maximum variation of the perturbation amplitude |U | in the lowest approximation can be expressed as: max ||U/U0| − 1| < 2 ∣∣d2 − dK2 + 8α0∣∣ /Ω2, (6.64) and by increasing Ω can be made arbitrarily small. In all of this the crucial point 144 Figure 6.8: Maximum amplitude of oscillations plotted against Ω andK. Parameters used: β0 = 0, β1 = 1, d = 8/3, (a) α0 = 0.5, a0 = 0, (b) α0 = 0.05, a0 = 0 and (c) α0 = 0.05, a0 = 0.15. The contour boundaries, starting from the color purple, are: 2, 5, 10, 20, 50, 100. is that the amplitude, while oscillating, remains finite. Hence, no collapse of the solitons in this case occurs. The solitons are modulationally stable against the perturbations considered. 145 Chapter 7 Conclusion Here we will give a brief summary of the results obtained in this Thesis. In Chapter 2, using a new ansatz developed in (143) we construct to the best of our knowledge the first exact spatiotemporal travelling wave and soliton solutions for the (2+1)D and (3+1)D nonlinear Schro¨dinger equations. For the first time, the role of chirp has been clearly displayed for multidimensional NLSEs. We obtain solutions to the (3+1)D NLSE both with and without chirp. The Jacobi elliptic function allows a continuous transformation of periodic solutions into solitary solution with the variation of parameter M . We then generalized our method to obtain for the first time solutions for the case of normal dispersion. The method is generalized further to cover nonlinearities of arbitrarily high order, although additional limitations are placed on our solutions so only for a handful of very specific parameters are solutions obtained. In Chapter 3, we use the ansatz developed in Chapter 2 to find exact solutions of the Gross-Pitaevskii equation. One discovers that finding exact solutions using our method is connected to finding the solution of the Ricatti equation. Assuming constant strength of field and diffraction/dispersion coefficient, one obtains solutions that either decay or blow up, remaining stable only for a critical value of gain. However, for the strength of quadratic field and diffraction/dispersion coefficient of sinusoidal form one obtains stable spatiotemporal travelling wave and soliton solutions. 146 Additional solutions are also found through a detailed exploration of the Ricatti equation. A general method of obtaining a wide class of solutions to the Ricatti equation is developed. These results are used to construct solutions for several new systems including those with Feschbach resonance. In Chapter 4, we obtain solutions for the related system of the NLSE within a linear potential. One obtains various solutions including ones that combine two periodic sine functions when the strength of linear field and diffraction/dispersion coefficient are of sinusoidal form. In Chapter 5, we apply our ansatz in a generalized form to the case of two co- or counter-propagating beams. We obtain novel solutions for the case where the coupling constant is equal to c = 3 both with and without chirp. The use of the Jacobi dn function produces solutions which don’t have any singularities. Finally, in Chapter 6 we analyze the stability of solutions given in Chapter 2 for the (3+1)D NLSE with normal or anomalous dispersion and the (2+1)D time- independent NLSE. For the (2+1)D solutions, we obtain stability for dark solitons and conditional stability for bright solitons, meaning we need to apply dispersion management to keep the solitons stable, i.e. the diffraction/dispersion coefficient must oscillate around 0. The management function dynamically stabilizes the non- linear structure of the transversal perturbation of the soliton if its mean value is zero. Reducing the period of the management function can be achieved by arbitrary limitations of the perturbation level. If, however, the mean value is different from zero, the amplitude modulation perturbation exponentially increases in the case of a solution without chirp and linearly in the case with chirp. In both cases we have the Lyapunov instability. For the (3+1)D case we obtain stability for the temporal bright solitons for normal dispersion and dark solitons for anomalous dispersion. All other types of solitons are conditionally stable. For the Gross-Pitaevskii equation we obtain that dark solitons are always stable and bright solitons are always conditionally stable. The obtained results are verified using computer simulations. The results obtained in this Thesis are in general amendable to further general- izations and applications to novel systems. Of particular interest is to generalize the 147 techniques used in this Thesis to other functions, for example the Weierstrass elliptic function. The stability of the solutions establishes the possibility of experimental verification sometime in the future. 148 Bibliography [1] R. W. Boyd, ”Nonlinear Optics”, Boston: Academic Press (1992). [2] M. A. Karpierz and G. I. Stegeman, Photonics Letters of Poland, 1(4), 145 (2009). [3] N. N. Akhmediev and A. Ankiewicz, ”Solitons: Nonlinear pulses and beams”, Chapman and Hall, London (1997). [4] Y. S. Kivshar and G. P. Agrawal, ”Optical solitons, from fibers to photonic crystals”, Academic Press, New York (2003). [5] C. S. Gardner et al., Phys. Rev. Lett. 19, 1095 (1967). [6] V. E. Zakharov and A. B. Shabat, Sov. Phys. JETP. 34, 62 (1972). [7] P. L. Sadchev, ”Self-similarity and beyond”, Chapman & Hall/CRC (2000). [8] R. Hirota, J. Math. Phys. 14, 810 (1973). [9] J. Satsuma, J. Phys. Soc. Jpn. 46, 359 (1979). [10] M. P. Miura, ”Backlund Transformation”, Springer-Verlag, Berlin (1978). [11] F. Cariello and M. Tabor, Physica D 53, 59 (1991). [12] Z. Fu et al., Phys. Lett. A 290, 72 (2001). [13] J. Kerr, Phil. Mag. 50, 337-348 and 446-458 (1875). [14] F. Pockels, Mathematisch-Physikalische Klasse 39, 1 (1894). [15] M. Go¨ppert-Mayer, Ann. Phys. 9, 273 (1931). [16] A. D. Buckingham, Proc. Phys. Soc. B 69, 344 (1956). 149 [17] A. Piekara and S. Kielich, J. Chem. Phys. 29, 1297 (1958). [18] T. H. Maiman, Nature 187, 493 (1960). [19] P. A. Franken et al., Phys. Rev. Lett. 7, 118 (1961). [20] R. W. Terhune, P. D. Maker, and C. M. Savage, Phys. Rev. Lett. 8, 404 (1962). [21] M. Bass, et al., Phys. Rev. Lett. 8, 18 (1962). [22] M. Bass et al., Phys. Rev. Lett. 9, 446 (1962). [23] W. Kaiser and C. G. B. Garett, Phys Rev. Lett. 7, 229 (1961). [24] I. D. Abella, Phys. Rev. Lett. 9, 453 (1962). [25] P. D. Maker, R. W. Terhune and P. C. Savage, Phys. Rev. Lett. 12, 507 (1964). [26] P. D. Maker and R. W. Terhune, Phys. Rev 137 (3 A), A801 (1965). [27] E. J. Woodbury and W. K. Ng, Proc. IRE 50, 2347 (1962). [28] G. Mayer and R. Gires, C. R. Acad. Sci. 258, 2039 (1963). [29] R. R. Alfano and S. L. Shapiro, Phys. Rev. Lett. 24, 592 (1970). [30] R. H. Stolen and C. Lin, Phys. Rev. A 14, 1448 (1978). [31] M. N. Islam, et al., Opt. Lett. 12, 625 (1987). [32] R. L. Carman et al., Phys. Rev. Lett. 17 (26), 1281 (1966). [33] R. H. Stolen, IEEE J. Quantum Electron. 11, 100 (1975). [34] Z. Y. Ou et al., Opt. Lett. 17, 640 (1992). [35] J.-P. Meyn and M. M. Fejer, Opt. Lett., Vol. 22, 1214 (1997). [36] A. Yariv, IEEE J. Quantum Electron. QE-14, 650 (1978). [37] G. Unnikrishnan, J. Joseph, and K. Singh, Appl. Opt. 37, 8181 (1998). [38] L. Yang et al., Commun. Theor. Phys. 49, 107 (2008). [39] K. C. Kao and G. A. Hockham, Proc. IEE 113, 1151 (1966). 150 [40] F. P. Kapron, D. B. Keck and R. D. Maurer, Appl. Phys. Lett. 17, 423 (1970). [41] T. Miya et al., Electron. Lett. 15, 106 (1979). [42] G. P. Agrawal, ”Fiber-Optic Communication Systems”, 2nd ed., Wiley, New York (1997). [43] D. R. Smith and N. Kroll, Phys. Rev. Lett. 85, 2933 (2000). [44] D. R. Smith, J. B. Pendry and M. C. K. Wiltshire, Science 305, 788 (2004). [45] V. G. Veselago, Sov. Phys. Usp. 10, 509 (1968). [46] A. A. Zharov, I. Shadrivov and Y. Kivshar, Phys. Rev. Lett. 91, 037401 (2003). [47] S. Tang et al., Opt. Exp. 19, 18283 (2011). [48] C. Sulem and P. Sulem, ”The Nonlinear Schro¨dinger Equation: Self-Focusing and Wave Collapse”, Springer-Verlag, Berlin (2000). [49] L. Gagnon and P. Winternitz P., J. Phys. A 21, 1493 (1988). [50] L. Gagnon and P. Winternitz P., J. Phys. A 22, 469 (1989). [51] F. O¨. Idlay, ”Derivation of the General Propagation Equation”, Retrieved from http://www.fen.bilkent.edu.tr/∼ idlay/courses/2006/577/NLSE deriva- tion.pdf [52] O. Cohen et al., Opt. Lett. 31, 954 (2006). [53] A. D. Buckingham and B. J. Orr, Rev. Chem. Soc. 21, 195 (1967). [54] V. Manakov, Sov. Phys. JETP 38, 248 (1974). [55] N. N. Akhmediev et al., Phys. Rev. Lett. 81, 4632 (1998). [56] J. S. Russell, ”14th meeting of the British Association Reports”, 311-390, York (1844). [57] E. Fermi, J. Pasta, and S. Ulam, Studies of Nonlinear Problems, Los Alamos National Laboratory, Technical Report, Document LA-1940 (1955). [58] N. J. Zabusky and M. D. Kruskal, Phys. Rev. Lett. 15, 240 (1965). 151 [59] A. Hasegawa and F. Tappert, Appl. Phys. Lett. 23, 142 (1973). [60] A. Hasegawa and F. Tappert, Appl. Phys. Lett. 23, 171 (1973). [61] D.-S. Wang, ”Optical Solitons in a Nonlinear Fiber Medium with Higher-Order Effects, Recent Progress in Optical Fiber Research”, M. Yasin (Ed.), ISBN: 978- 953-307-823-6, InTech. Retrieved from: www.intechopen.com/books/recent- progress-in-optical-fiber-research/optical-solitons-in-a-nonlinear-fiber- medium-with-higher-order-effects (2012). [62] L. F. Mollenauer, R. H. Stolen, and J. P. Gordon, Phys. Rev. Lett. 45, 1095 (1980). [63] P. Emplit, M. Haelterman, and J.-P. Hamaide, Opt. Lett. 18, 1047 (1993). [64] A. Hasegawa, ”Optical Solitons in Fibers”, Springer-Verlag, Berlin (1989). [65] J. E. Bjorkholm and A. Ashkin, Phys. Rev. Lett. 32, 129 (1974). [66] A. Barthelemy, S. Maneuf and C. Froehly, Opt. Commun. 55, 201 (1985). [67] J. S. Aitchison et al., Electron. Lett. 28, 1879 (1992). [68] J. Beeckman et al., Opt. Express 12, 1011 (2004). [69] G. Assanto, M. Peccianti and C. Conti, Opt. Photon. News 14, 44 (2003). [70] U. A. Laudyn et al., Photonics Lett. of Poland 1, 7 (2009). [71] A. I. Strinic et al., Opt. Express 14, 12310 (2006). [72] Y. V. Izdebskaya et al., Opt. Express 18, 3258 (2010). [73] D. N. Christodoulides and R. I. Joseph, Opt. Lett. 13, 53 (1988). [74] Y. S. Kivshar, Opt. Lett. 17, 1322 (1992). [75] S. T. Cundiff et al., Phys. Rev. Lett. 82, 3988 (1999). [76] D. Y. Tang et al., Phys. Rev. Lett. 101, 153904 (2008). [77] A. S. Desyatnikov, Y. Kivshar and L. Tomer, Progress in Optics 47, 291 (2005). 152 [78] Y. Kivshar and D. Pelinkovsky, Phys. Rep. 331, 117 (2000). [79] C.-C. Jeng et al., Phys. Rev. Lett. 92, 043904 (2004). [80] I. Towers and B. A. Malomed, J. Opt. Soc. Am. B 19, 537 (2002). [81] I. Yakimenko et al., Phys. Rev. E 73, 066605 (2006). [82] K. Motzek et al., J. Opt. Soc. Am. B 22, 1437 (2005). [83] C. A. Swartzlander Jr. and C. T. Law, Phys. Rev. Lett. 69, 2503 (1992). [84] V. Tikhonenko, V. Lashkin and O. Prikhodko, Opt. Lett. 21, 1129 (1996). [85] J. C. Neu, Physica D 43, 385 (1990). [86] A. Desyatnikov, A. Sukhorukov and Y. Kivshar, Phys. Rev. Lett. 95, 203904 (2005). [87] S. Lopez-Aguayo, A. Desyatnikov and Y. Kivshar, Opt. Exp. 14, 7903 (2006). [88] Y. Zhang et al., Opt. Exp. 18, 27846 (2010). [89] J. W. Fleischer et al., Phys. Rev. Lett. 92, 123904 (2004). [90] B. Terhalle et al., Phys. Rev. A 79, 043821 (2009). [91] M. Belic´ et al., Phys. Rev. E 65, 066610 (2002). [92] J. J. Garcia-Ripoll et al., Phys. Rev. Lett. 85, 82 (2000). [93] K. Motzek et al., Phys. Rev. E 68, 066611 (2003). [94] M. Belic´ et al., Opt. Exp. 12, 708 (2004). [95] M. Belic´ et al., Opt. Exp. 14, 794 (2006). [96] D. Jovic´ et al., Opt. Exp. 13, 4379 (2005). [97] T. Aktosun, F. Demontis and C. van der Mee, Inverse Problems 23, 2171 (2007). [98] P. Drazin and R. Johnson, ”Solitons: An Introduction”, Cambridge University Press (1989). 153 [99] M. Lax, J. H. Butteh and G. P. Agrawal, J. Appl. Phys. 52, 109 (1981). [100] V. E. Zakharov and S. V. Manakov, Journal of Theoretical and Mathematical Physics 19, 551 (1974). [101] M. J. Ablowitz et al., Stud. Appl. Math. 53, 249 (1974). [102] L. Debnath, ”Nonlinear Partial Differential Equations”, Third Edition, Birkha¨user (2010). [103] T. Aktosun, arXiv:0905.4746 (2009). [104] T. Busse, Ph.D. thesis, University of Texas at Arlington (2008). [105] M. Soljacˇic´, M. Segeev and C. R. Menzuk, Phys. Rev. E 61, 1048 (2000). [106] Y. Ye, Mathematical Problems in Engineering, 298980 (2009). [107] B. Yang, W.-P. Zhong and M. Belic´, Commun. Theor. Phys. 53, 937 (2010). [108] S.-L. Xu, J.-C. Liang and Z.-M. Li, Chin. Phys. B 20, 114214 (2011). [109] Y. Xia, W.-P. Zhong and M. Belic´, Acta Physica Polonica B 42, 1881 (2011). [110] H.-P. Zhu, Commun. Theor. Phys. 58, 67 (2012). [111] N. Z. Petrovic´, M. Belic´ and W.-P. Zhong, Phys. Rev. E 83, 026604 (2011). [112] K. Senthilnathan, A. M. Abobaker and K. Nakkeeran, PIERS Proceedings, Moscow, Russia, 1396 (2009). [113] P. P. Goldstein, Acta Physica Polonica A 112, 1071 (2007). [114] Z.-D. Li et al., Annals of Physics 322, 2545 (2007). [115] Q.-Y. Li et al., Opt. Comm. 282, 1676 (2009). [116] R. Hirota, ”The Direct Methods in Soliton Theory”, Cambridge University Press, Cambridge (2004). [117] W. T. Wu, Sci. Sinica 21, 159 (1978). [118] W. T. Wu, Journal of Sys. Sci. Math. Scis. 4, 207 (1984). 154 [119] L. Yang, J. Liu, and K. Yang, Phys. Lett. A 278, 267 (2001). [120] Z. Fu, S. Liu and S. Liu, Phys. Lett. A 299, 507 (2002). [121] Z. Yan, J. Phys. A 36, 1961 (2003). [122] M. Inc and D. J. Evans, International Journal of Computer Mathematics 81, 313 (2004). [123] Z. Yan and H.Q. Zhang, Phys. Lett. A 285, 355 (2001). [124] Y. B. Zhou, M. L. Wang, and Y. M. Wang, Phys. Lett. A 308, 31 (2003). [125] X.-Y. Li, X.-Z. Li and M.-L. Wang, Commun. Theor. Phys. 45, 9 (2006). [126] B. Hong and D. Lu, International Journal of Nonlinear Science 7, 360 (2009). [127] Y. Chen and Q. Wang, Chaos, Solitons & Fractals 24, 745 (2005). [128] Y. Chen, Q. Wang and B. Li, Z. Naturforch. 59a, 529 (2004). [129] Z. Yan, Computer Physics Communications 153, 145 (2003). [130] A. T. Ali, Journal of Computational and Applied Math 235, 4117 (2011). [131] B. Hong, D. Lu and F. Sun, World Journal of Modeling and Simulation 5, 216 (2009). [132] E. J. Parkes, B. R. Duffy and P. C. Abbott, Phys. Lett. A 295, 280 (2006). [133] Y.-J. Ren and H.-Q. Zhang, Chaos, Solitons & Fractals 27, 959 (2006). [134] S. Zekovic´ et al., arXiv: 1212.0374 (2012). [135] S. Zhang, T.-C. Xia, Applied Mathematics & Computation 183, 1190 (2006). [136] Y. B. Zhou, M. L. Wang, and T. D. Miao, Phys. Lett. A 323, 77 (2004). [137] A. Ebaid and S. M. Khaleed, Journal of Computational and Applied Mathe- matics 235, 1984 (2011). [138] J.-M. Zhu et al., Chin. Phys. 13, 798 (2004). [139] S. N. Alsaedi, Int. J. Contemp. Math. Sci. 5, 341 (2010). 155 [140] V. I. Kruglov, A. C. Peacock, and J. D. Harvey, Phys. Rev. E 71, 056619 (2005). [141] M. Abramovitz and I. A. Stegun, eds., ”Handbook of Differential Equations”, Chapter 16, Dover (1965). [142] M. Belic´, N. Z. Petrovic´, W.-P. Zhong, R. H. Xie and G. Chen, Phys. Rev. Lett. 101, 123904 (2008). [143] W.-P. Zhong, R.-H. Xie, M. Belic´, N. Z. Petrovic´, G. Chen and L. Yi, Phys. Rev. A 78, 023821 (2008). [144] N. Z. Petrovic´, M. Belic´, W.-P. Zhong, R.-H. Xie and G. Chen, Opt. Lett. 34, 1609 (2009). [145] N. Z. Petrovic´, M. Belic´ and W.-P. Zhong, Phys. Rev. E 81, 016610 (2010). [146] N. Z. Petrovic´, N. Aleksic´, A. Al Bastami and M. Belic´, Phys. Rev. E 83, 036609 (2011). [147] A. Al Bastami, N. Z. Petrovic´ and M. R. Belic´, Electron. J. Diff. Eqs., Vol. 2010, No. 66, 1 (2010). [148] A. Al Bastami, M. R. Belic´, D. Milovic´ and N. Z. Petrovic´, Phys. Rev. E 84, 016606 (2011). [149] N. Z. Petrovic´, H. Zahreddine and M. Belic´, Phys. Scr. 83, 065001 (2011). [150] N. Z. Petrovic´ and H. Zahreddine, Phys. Scr. T149, 014039 (2012). [151] N. N. Akhmediev and A. B. Shabat, ”Solitons: Nonlinear pulses and wavepack- ets”, Chapman and Hall, London, (1997). [152] L. Berge´, et al., Opt. Lett. 25, 1037 (2000). [153] B. A. Malomed, ”Soliton Management in Periodic Systems”, Springer, New York (2006). [154] G. D. Montesinos, et al., Chaos 15, 033501 (2005). [155] H. Saito and M. Ueda, Phys. Rev. Lett. 90, 040403 (2003). 156 [156] W.-P. Zhong and L. Yi, Phys. Rev. A 75, 061801 (2007). [157] B. A. Malomed, et al., J. Opt. B 7, R53 (2005). [158] M. Gedalin, T. C. Scott, and Y. B. Band, Phys. Rev. Lett. 78, 448 (1997). [159] S. Konar and A. Biswas, Opt. Quant. Electron. 36, 905 (2004). [160] F. Dalfovo et al., Rev. Mod. Phys. 71, 463 (1999). [161] L. P. Pitaevskii, and S. Stringari, ”Bose-Einstein Condensation”, Oxford Uni- versity Press, Oxford (2003). [162] E. P. Gross, Phys. Rev. 106, 161 (1957). [163] E. P. Gross, Il Nuovo Cimento 20, 454 (1961). [164] V. L. Ginzburg and L. P. Pitaevskii, Sov. Phys. JETP 7, 858 (1958). [165] L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 (1961). [166] N. Parker, University of Durham, Ph. D. Thesis. [167] A. Hasegava and M. Matsumoto, ”Optical Solitons in Fibers”, Springer (2003). [168] W. Bao, D. Jaksch, and P. A. Markowich, Journal of Computational Physics 187, 318 (2003). [169] S. Adhikari and P. Muruganandam, J. Phys. B: At. Mol. Opt. Phys. 35, 2831 (2002). [170] S. Burger et al., Phys. Rev. Lett 83, 5198 (1999). [171] H. Sakaguchi and B. A. Malomed, Phys. Rev. E 73, 026601 (2006). [172] R. Atre, P.K. Panigrahi, and G.S. Agarwal, Phys. Rev. E 73, 056611 (2006). [173] Q. Yang and H.-J. Zhang, Chinese J. Phys. 46, 457 (2008). [174] S. K. Adhikari, Phys. Rev. A 69, 063613 (2004). [175] S. K. Adhikari, Phys. Rev. E 71, 016611 (2005). 157 [176] J. F. Riccati, ”Actorum Eruditorum quae Lipsiae Publicantur”, Suppl. 8, 66 (1724). [177] E. L. Ince, ”Ordinary Differential Equations”, Dover Publications, New York (1956). [178] H. Davis, ”Introduction to nonlinear differential and integral equations”, Courier Dover Publications (1962). [179] W. T. Reid, ”Riccati Differential Equations”, Academic Press, New York (1972). [180] E. Kamke, ”Differentialgleichungen: Lo¨sungsmethoden und Lo¨sungen, I, Gewo¨hnliche Differentialgleichungen”, B. G. Teubner, Leipzig (1977). [181] D. Zwillinger, ”Handbook of Differential Equations”, Academic Press, Boston (1989). [182] A. D. Polyanin and V. F. Zaitsev, ”Handbook of Exact Solutions for Ordi- nary Differential Equations”, 2nd Edition, Chapman & Hall/CRC, Boca Raton (2003). [183] F. Schwabl, ”Quantum Mechanics”, Springer, Berlin (1992). [184] E. S. Fraga, ”The Schro¨dinger and Riccati Equations”, Lecture Notes in Chem- istry Vol. 70, Springer, Berlin (1999). [185] R. Shankar, ”Principles of Quantum Mechanics”, Plenum, New York (1980). [186] R. Liboff, ”Introductory Quantum Mechanics”, 4th ed., Addison Wesley, New York (2002). [187] A. Khare and U. Sukhatme, ”Supersymmetry in Quantum Mechanics”, World Scientific, Singapore (2001). [188] C. Yuce and A. Kilic, Phys. Rev. A 74, 033609 (2006). [189] S. Inouye, Nature 392, 151 (1998). [190] H.-H. Chen and C.-S. Liu, Phys. Rev. Lett. 37, 693 (1976). 158 [191] Q. Yang, J.-F. Zhang, Opt. Comm. 258, 35 (2006). [192] M. Petrovic´ et al., Laser & Photonics Reviews 5, 214 (2011). [193] G. P. Agrawal, Phys. Rev. Lett. 59, 880 (1987). [194] Z. Rapti et al., Phys. Rev. E 69, 017601 (2004). [195] J. A. C. Weideman and B. M. Herbst, Siam. J. Numer. Anal. 23, 485 (1986). [196] M. Lax et al., J. Opt. Soc. Am. A 2, 731 (1985). [197] J. Javanainen and J. Ruostekoski, J. Phys. Mat. Gen. 39, L179 (2006). [198] G. M. Muslu and H. A. Erbay, Mathematics and Computers in Simulation 67, 581 (2005). [199] J. W. Cooley and J. W. Turkey, Math. Comput. 19 (90), 297 (1965). [200] Y. Castin, ”Coherent atomic matter waves”, R. Kaiser, C. Westbrook, and F. David, eds., Les Houches Session LXXII, Springer (2001). 159 160 Биографија Никола Петровић је рођен 12. 03. 1980. године у Београду. Завршио је Математичку Гимназију 1999. године као ученик генерације са просеком 5.00. Дипломирао је физику и математику у јуну 2003. године на Масачусетс институту за технологију (Massachusetts Institute of Technology) са просеком 4.5 (на скали од 0 до 5). Дипломски рад је био на тему кодова за исправљање грешака у квантним компјутерима: “Constructing an Infinite Class of Perfect Codes”, са оценом Б (9). Ментор му је био проф. Исак Чуанг (Isaac Chuang). Објавио је са још три коаутора књигу “The IMO Compendium”, са свим задацима предложеним на Међународним математичким олимпијадама (Спрингер-Верлаг, Берлин, 2006). Пријавио се 2004. године на Београдски Универзитет на магистарске студије, а 2008. прешао са магистарских на докторске студије. Положио је све испите из изборних предмета на докторским студијама са просечном оценом 10. На седници Научно-наставног већа одржаној 7.12.2009. године је одобрена израда његове докторске дисертације. Од 2004. године Никола Петровић је у радном односу са Институтом за Физику у Београду. Његов статус је замрзнут од августа 2005. године када одлази на Тексашки А&М универзитет у Катару (Теxаs А&М University at Qatar) где је био запослен као лабораторијски координатор и радио такође као асистент све до јула 2012. године када се враћа у Институт за Физику. У септембру 2012. је изабран у звање истраживача сарадника. До треунтка одбране дисертације публиковао је 13 радова у међународним часописима. Ожењен је супругом Ташаном Петровић и има сина Бориса. 161 Biography Nikola Petrović was born on 12.03.1980 in Belgrade. He finished the High School of Mathematics as a valedictorian with an average of 5 (A). He graduated in Physics and in Mathematics in June 2003. at the Massachusetts Institute of Technology with an average of 4.5 (on a scale from 0 to 5). The topic of his undergraduate Thesis was quantum error correction: “Constructing an Infinite Class of Perfect Codes” and he received a grade of B. His mentor was prof. Isaac Chuang. He has published with three other co-authors the book “The IMO Compendium”, a collection of all problems shortlisted for the International Mathematical Olympiad (Springer-Verlag, Berlin 2006). He applied in 2004. for a Master’s Degree at Belgrade University, switching over to pursuing a Ph. D. in 2008. He has passed all post-graduate classed with an average grade of 10 (A). On the faculty meeting on 7.12.2009 the topic for his Thesis was approved. Since 2004. He is employed at the Institute of Physics in Belgrade. He was placed on leave on August 2005. when he left for Теxаs А&М University at Qatar) to work as a laboratory coordinator and a teacher associate (TA) until June 2012. when he returned to the Institute of Physics in Belgrade. In 2012. he became a research associate of the institute. Prior to the defence of his Thesis he has published 13 papers in international journals. He is the married to Tašana Petrović and has a son Boris.