KoKoJevy Numerické řešení Grossovy-Pitaevského rovnice David Linhart Naším úkolem je numerické řešení Grossovy-Pitaevského rovnice, která přibližně popisuje kondenzát slabě interagujících atomů v harmonické pasti. Pro účely numerického řešení si ji zapíšeme v bezrozměrných veličinách a navíc budeme určovat pouze základní stav, který je sféricky symetrický:1 − d2 dξ2 + ξ2 + 8πA|ϕ(ξ)|2 u(ξ) = ˜µu(ξ), (1) kde ξ = r/aHO, A = N a aHO je parametr charakterizující poměr odpudivé interakční a kinetické energie, ˜ϕ(ξ) = ϕ√ Na3 HO = u(ξ) ξ , ˜µ = µ 1/2ℏωHO a aHO = ℏ mωHO je charakteristická délka. Jedná se o numerické řešení okrajové úlohy, vzhledem k tomu, že operátor závisí na hustotě atomů a hustota atomů přirozeně závisí na jejich vlnové funkci, jedná se o selfkonzistentní problém, nemůžeme tudíž použit připravené rutiny pro řešení okrajových problémů. V našem řešení použijeme diskrétní množinu hodnot ξ, derivaci v rovnici převedeme na diferenci a následně iterativně dopočítáváme hodnotu potenciálu V = 8πA|˜ϕ(ξ)|2. Pro převedení derivace na diferenci využijeme aproximativního vztahu: − d2u dξ2 ≈ (2uj − uj+1 − uj−1) 1 (ξmax/n)2 . Využitím tohoto vztahu se problém změní na problém diagonalizace tridiagonální matice, která má na diagonále hodnoty: 2 (ξmax/n)2 + ξ2 + A|ϕ(ξ)|2 a pod a nad diagonálou −1 (ξmax/n)2 . Vlnová funkce je poté obsažena ve vektoru příslušícímu nejnižší vlastní hodnotě. Při každé iteraci je nutné normovat vlnovou funkci, aby stále platilo, že: 1 = ˜ϕ2 (ξ)d3 ξ = u2 ξ2 4πξ2 dξ ≈ u2 4π ξmax n . (2) Pro velké hodnoty A je nutné stabilizovat iterační postup tlumením změn potenciálu podle předpisu: V (n+1) = (1 − λ)V (n) + λA(˜ϕ(n))2. Z výsledného grafu je zřejmé, že jak se zvětšuje odpudivá interakční energie oproti kinetické energii, tak je obláček kondenzátu rozprostřenější v prostoru. Zároveň roste i hodnota chemického potenciálu - bude těžší atomy přidat do kondenzátu, tedy silnější odpudivá interakce způsobí větší energii systému, jak jsme obecně zvyklí. 1 Detaily úpravy jsou uvedené v zadání. Ukázka řešení v jazyce Julia pomocí balíčku LinearAlgebra § ¤ 1 using LinearAlgebra 2 3 n=100 4 ξmax=6 5 ξ=Array(LinRange(6/n,ξmax,n)) 6 DiffFact=(ξmax/n)ˆ2 7 NormFact=4*π*ξmax/n 8 λ=1/12 9 dlu=-ones(n-1)/DiffFact 10 11 Alist=[0,.1,1,10,100] 12 for A ∈ Alist 13 u=ones(n) 14 ϕ=u./ξ 15 ϕ/=sqrt(sum(u.ˆ2 .*NormFact)) 16 V=zeros(n) 17 for i=1:250 18 ϕlast=ϕ 19 20 DiaEle = (ξ.ˆ2 .+ V .+ (2/DiffFact)) 21 operator=Matrix(Tridiagonal(dlu,DiaEle,dlu)) 22 eig=eigen(operator) 23 u=eig.vectors[:,1] 24 global µ=eig.values[1] 25 26 ϕ=u./ξ 27 ϕ/=sqrt(sum(u.ˆ2 .*NormFact)) 28 V=(1-λ).*V+λ*A .*abs2.(ϕ)*8*π 29 end 30 end ¦ ¥ 2