TVE-F 19002
Examensarbete 15 hp Juni 2019
NON-THERMAL MODELING OF ENERGY PROPAGATION
CARRIED BY PHONON AND MAGNONS David Dahlgren
Amela Mehic
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Abstract
Non-thermal modeling of energy propagation carried by phonons and magnons
David Dahlgren, Amela Mehic
Heat transport by phonons and magnons in crystals due to a local perturbation intemperature is described by the Boltzmann transport equation. In this project a onephonon one magnon system was studied in a one dimensional rod with reflectiveboundaries. Using the splitting algorithm the problem was reduced to a transport andcollision term. The resulting stabilization time for a initial phonon and magnondistribution and the respective temperatures at equilibrium was calculated.
This studyshows how energy propagates in crysials and gives further understanding of how thecoupling of phonons and magnons affect heat transport.
Ämnesgranskare: Chao Xu
Handledare: Pablo Maldonado, Jerome Hurst
1 Popul¨ arvetenskaplig sammanfattning
I v˚art samh¨alle anv¨ands elektronik i stor utstr¨ackning och utan den skulle inte v˚art samh¨alle fungera. Elektronikens utveckling ¨ar dels beroende av komponen- ternas materialegenskaper. En sv˚arighet ¨ar upphettningen av komponenterna under l¨angre perioder av anv¨andning som har visat sig p˚averka deras egenskaper och i vissa fall verkningsgrad. D¨arf¨or ¨ar det extremt relevant att studera hur metaller och kristaller p˚averkas av extrema temperaturf¨or¨andringar.
I detta projekt har v¨armeledning i en endimensionell kristallfilm som uts¨atts f¨or en lokal temperaturf¨or¨andring studerats. V¨armeledning p˚a mikroskopisk niv˚a beskrivs av fononer- och magnonkopplingars propagering och interaktion med varandra. Fononer och magnoner ¨ar quasi partiklar som beskriver vibra- tioner och spin i en gitterstruktur med hj¨alp av en transformation till ett annat rum som g¨or det m¨ojligt att diagonalisera Hamiltonianen. Genom att anv¨anda Boltzmanns transport ekvation f¨or att beskriva hur dessa fononer och magnoner transporteras och interagerar med varandra kan v¨armeledning i kristallen beskri- vas. I v˚art projekt studerades en fonon en magnon system.
F¨or att numeriskt l¨osa Boltzmann ekvationen anv¨andes en ”splitting algo- rithm” som delar upp problemet i en transport och en collision term. Genom att simulera upprepande ¨over olika l˚anga tidssteg noterades tiden det tog f¨or sys- temet att stabiliseras. Efter stabiliseringstiden hade fononerna och magnonerna en respektive konstant homogen temperatur som ber¨aknades. Slutreslutatet visade att det teoretiska och simulerade v¨ardena st¨amde ¨overens och gav oss en inblick i hur en extrem, lokal temperaturf¨or¨andring p˚averkar stabiliseringstiden hos en kristall.
Contents
1 Popul¨arvetenskaplig sammanfattning 3
2 Introduction 5
2.1 Goals . . . . 5
2.2 Structure . . . . 5
3 Theory 6 3.1 The Classical Problem and Dispersion Relation . . . . 6
3.2 Phonons . . . . 8
3.2.1 Hamiltonian For a Chain Of Atoms . . . . 8
3.2.2 Using Fourier Transform to Simplify the Hamiltonian . . 8
3.2.3 Express the Hamiltonian as a Harmonic Oscillator Using Creation And Annihilation Operators . . . . 9
3.3 Magnons. . . . 10
3.3.1 Heisenberg Hamiltonian and Dispersion Relation For Magnons 10 3.4 The Coupling Of Phonon and Magnons . . . . 12
3.5 The Boltzmann Equation . . . . 12
3.5.1 Equilibrium Temperature for Phonons and Magnons . . . 14
3.6 Boundary and initial conditions for solving the Boltzmann equation 14 4 Computational Method 15 4.1 Splitting Algorithm. . . . 15
4.1.1 The Transport Term . . . . 16
4.1.2 The Collision Term. . . . 16
4.1.3 Ratio of Minimum and Maximum Values . . . . 17
5 Results 17 5.1 Numerical Values For Constants . . . . 17
5.2 Simulation With Initial Excited Phonon and Magnon. . . . 18
6 Discussion 22 7 Conclusion 23 7.1 Appendix . . . . 25
2 Introduction
How different forms of energy, such as heat, is transported in a crystal structure has become a central question in condensed matter science, due to the strong out-of-equilibrium conditions and different scatterings involved in the process.
Additionally, to gain a better understanding of this process is essential for de- velopment of modern technologies. In this project we will address the question by using the Boltzmann transport equation considering not only the out-of- equilibrium conditions but also the coupling between phonons and magnons.
To simulate such extreme conditions, in our work we will use a spike in tem- perature between two point on a small one dimensional, 4nm long, rod. This temperature excitation will supply energy to the system and excite phonon and magnons. Then we will use the Boltzmann equation to describe how the phonon- magnon interaction affects the energy dispersion. The problem will be solved numerically by implementing a code in C++. Utilizing the splitting algorithm method the complicated problem will be reduced to a transport and collision differential equation. The input parameters necessary to solve the equation are group velocities, decay constant, and energies for one phonon and one magnon mode. The system dynamics will be simulated and will provide answers to what happens at different times. We will also investigate the time needed for the sys- tem to stabilize for different initial conditions. Results of or simulations will be presented in plots from Matlab. Additionally numerical values at equilib- rium will be utilized to calculate the steady state temperature for phonon and magnons.
2.1 Goals
The goal of this project is to simulate phonon-magnon heat transport and dy- namics using the Boltzmann equation. Investigating the behaviour of magnon and phonon distributions throughout time and extracting the time needed for the system to stabilize. After the stabilization of the system the new equilibrium temperature is calculated for magnons and phonons respectively.
2.2 Structure
The first part of this paper is a theory section where phonons, magnons and the Boltzmann equation is explained in substantial detail. The second section describes the numerical method of simulating the phonon-magnon interaction.
The third section reports on the simulation of a one phonon one magnon system to observe the dynamics and energy propagation. In The final section we discuss the results and how the methods in this project could be improved to produce a more accurate simulation of heat transfer in crystals.
3 Theory
Heat transfer in insulating crystals is mainly carried by phonon and magnons.The
??-quantum description of the quasiparticles is given by the Hamiltonian of crys- tal lattice and will be generally derived in this project. The out-of-equilibrium transport will be described by the Boltzmann equation, for witch it is neces- sary to understand the scattering between quasiparticles and their transport.
Due to the complex nature of the interaction between phonons and magnons, this will be approximated to the problem of one phonon one magnon interac- tion. These quasi particles are derived using a Hamiltonian of a crystal lattice utilizing Fourier transform and diagonalize the potential[1],[2],[3].
3.1 The Classical Problem and Dispersion Relation
To understand the concept of phonons it is helpful to first examine the prob- lem from a classical point of view. Imagine a chain of atoms coupled together with springs and with a equlibrium distace ,a, between them, a linear lattice.
The springs will represent the Coulomb force acting on every atom when there is a deviation from equilibrium. Each atom have a known position, qs, and a momentum, ps[1]. To solve this problem one could utilize Hookes’s law which says that the force extorted on a point connected with a spring is linearly pro- portional to the deviation from equilibrium. The forces on an atom sited at s, can be described as the sum of forces extracted from each atom on the atom s, becoming
Fs=X
n
Cn(qs+n− qs) (1)
where Cn is the spring constant, n is a specific atom in the lattice and qs
represents the position of the mass point of interest. From Newtons second law this expression becomes
md2qs
dt2 =X
n
Cn(qs+n− qs). (2)
where m is the mass for each atom[1]. Equation 2 describes a correlation motion relation of the linear lattice chain. The equation indeed is continuous in time, but also depends on discrete lattice points in space, making it a discrete differential equation. If only small deviations from equilibrium is considered the wavelength of the lattice chain must be long. In comparison with the long wavelength the lattice points are tightly packed allowing an approximation of continuity. Equation2 can be extended
mddt2q2s = C1[(qs+1− qs) − (qs− qs−1)(aa)2] + C2[(qs+2− qs) − (qs− qs−2)(2a2a)2]+
C3[(qs+3− qs) − (qs− qs−3)(3a3a)2] + ...
(3)
where qs+n− qs can be approximated by dqs, na by dx and from translation symmetry Cn= C−n[1].Equation3can be simplified too
md2qs
dt2 =d2qs
dx2[C1
ma2+C2
m(2a)2+C3
m(3a)2+ ...].
Comparing this to the one dimensional wave equation md2q
dt2 = v2d2q
dx2 (4)
and with the approximations above the solution will be a plane wave that can be expressed like
qs= Ae−i(ωt−ksa). (5)
A is the amplitude of the wave, ω is the angular frequency, k =2πλ is the wave vector and a is the length between two sequent lattice points[1]. This indicates that the expected solution to equation 2 will have a wave solution as well.
Putting equation5into equation 1and simplifying with translation symmetry, ω can be expressed as a function of k
ω2(k) = 2 m
X
n>0
Cn(1 − cos(nka)). (6)
This is commonly known as the dispersion relation that describes the lattice vibrations depending on the wave vector[1]. One can see that the frequency maximum is bound meaning it can not increase infinitely when the wavelength decreases. Arguing that the coupling force decreases in strength when the in- teractions between two atoms increases, allow us to consider only nearest neigh- bours simplifying the problem to n = 1. Substituting this into equation6yields the nearest neighbour thus approximation:
ω =
√4C1 m
sin1 2ka
. (7)
Lattice point vibrations are atomic displacements that are described by wave functions. Since only the lattice points are considered in the wave function, everything between then is irrelevant for the description of the problem and can be ignored. This with equation6yields that the same frequency ω can be given for different wave vectors k0 = k +2πa . Since the wave vector and wavelength are inversely proportional one can state that different wavelengths can describe the same lattice vibration. Comparing these waves in a continuous spectrum it must be stated that they are different because of different wavelength. However by looking at only the lattice points their values are the same. All information can be drawn from this range which is called the first Brillouin zone. The second Brillouin zone is the repetition of the first and so on[1].
3.2 Phonons
In reality there are difficulties with the classical approach in solving this prob- lem. The uncertainty principle which states that one can not know the exact position and momentum at a given time. The probability of finding a particle in a region of space is different between the classical and quantum case. However solving problem using quantum mechanics the result will not only lead to a relatively simple expression for a complicated problem but will also lead to the definition of phonons, which are a quantization of vibrating modes in a lattice structure[1].
3.2.1 Hamiltonian For a Chain Of Atoms
Now the problem will be solved using the quantum physics approach. To express the Hamiltonian for the system of atoms, the first step is to write the potential energy. One atom will feel a potential from every other point in the system which will be proportional to the square of the distance between two atoms.
Utilizing this fact the potential for every atom can be expressed as a sum of potentials
V (qˆ 1, q2, ..., qN) =1 2
X
s
[1 2
X
s0
Css0(ˆqs− ˆqs0)2] (8) where all atom positions are operators, N are the number of particles and Css0 are the spring constants between atom s and s’[1]. The first summation runs over all atoms that influence one atom and the second summation runs over all calculated potentials from every atom. The factor 12 is needed because of symmetry, (ˆqs− ˆqs0)2 = (ˆqs0 − ˆqs)2 and to not have double counting in the energy. The Hamiltonian for this system can be expressed as
H =ˆ 1 2
X
s
(pˆ2s 2m+1
2 X
s0
Css0(ˆqs− ˆqs0)2) = 1 2
X
s
ˆhs. (9)
This Hamiltonian is a sum of coupled harmonic oscillators where each term,ˆhs, represents an individual Hamiltonian for each of the atoms. In order to con- tinue the derivation the function that describe the system needs to have peri- odic boundary conditions. Introducing this condition results in N independent modes existing in the fist Brillouin zone[1].
3.2.2 Using Fourier Transform to Simplify the Hamiltonian
The difficulty with the Hamiltonian ,equation9, is the cross terms that couples every atom to each other. However by changing from (ˆqs, ˆps) → ( ˆQk, ˆPk), by a Fouirer transform, the problem can be simplified[1]. Using the discrete Fourier transform gives the following expression for ˆqs, ˆps
ˆ qs=
r1 N
X
k
Qˆkeiksa (10)
ˆ ps=
r1 N
X
k
Pˆke−iksa (11)
which expresses the lattice vibrations as a summation over all possible fre- quencies in k space[1]. By also imposing a periodic boundary conditions, the atom Ns+ 1 is equivalent to the 1st atom, due to the lattice structure and the first Brillouin zone. This transformation removes the coupling term and makes it easier to describe the system. When applying equation10,11 to equation9 the Hamiltonian in k space becomes
H =ˆ 1 2m
X
k
Pˆ−kPˆk+1 2
X
k
mω2kQˆ−kQˆk (12) where ωk is the dispersion relation, equation 6[1]. The clear difference be- tween equation9and12is that the summation is only over one variable k where the previous Hamiltonian summed over s and s0. As a result of the transforma- tion, one can define the atomic vibrations as a linear combination of quantized collective excitations called phonons.
3.2.3 Express the Hamiltonian as a Harmonic Oscillator Using Cre- ation And Annihilation Operators
The final step is to express the Hamiltonian in terms of annihilation and cre- ation operators, ˆak and ˆa+k respectively. This can be obtained by defining the operators in the following form
ˆ
ak= B ˆQk+ iA ˆP−k, ˆa+k = B ˆQ−kiA ˆPk (13) where A = √ 1
2m~ωk and B = pmωk
2~ [1]. Utilizing the fact that equation 12 is invariant when changing k to −k, that the summation runs over all k inside the fist Brillouin zone and symetric around k = 0 gives
H =ˆ X
k
( ˆNk+1
2)~ωk (14)
where ˆNk= ˆa+kˆak[1]. The final expression looks just like the Hamiltonian for a harmonic oscillator and the different eigenenergies for a mode k with population nk is given by
Enk,k= (nk+1
2)~ωk. (15)
Hence by transforming from real space coordinates to reciprocal space coordi- nates and diagonalize the potential, the Hamiltonian results in a set of N uncou- pled harmonic oscillators with quantized energies. The quantum of vibrational energy is a phonon. In this case only an one-dimensional lattice was considered however for more complicated geometries the work path will be similar.
3.3 Magnons
The local magnetic moments of each atom, localized spins, couples in a similar manner as phonons. If the ground state of the magnetic moment, ˆS, is aligned along ˆz, ˆSz, an external perturbation of the spins will result in a deviation from equilibrium. The resulting spin projection onto ˆz will be different, less than the unperturbed spin along ˆz, and as a consequence rotational dynam- ics will be initiated along the z-axis, leading to different energies for different spin configurations. These dynamics can be modeled by using the Heisenberg Hamiltonian.
Figure 1: Illustrations of magnon out-of-equilibrium spin. The arrows spin around the circle representing alteration from equilibrium.
In figure 1, one can see the dynamics of description above. When spin is in equilibrium the spin is aligned along the ˆz-axis. If an energy perturbation is introduced the spin departs from the ˆz-axis as seen in the figure on the left respectively right.
3.3.1 Heisenberg Hamiltonian and Dispersion Relation For Magnons The Hamiltonian that describes the spin interaction in a ferromagnetic material is the Heisenberg Hamiltonian
H = Jˆ X
iδ
SˆiSˆi+δ− 2µBB0
X
i
Sˆiz (16)
where J is the ferromagnetic exchange constant, i is a summation over all the strides, δ represents all the nearest neighbors, µB is the Bohr magneton and B0 is an external magnetic field along the ˆz-axis. This describe the exited interaction explained above and when all the spins are aligned in the ˆz direction the system is in equilibrium. It is also important to note that
[ ˆSx, ˆSy] = i ˆSz (17) and therefore we could have chosen ˆx and ˆx as reference axis[3]. The next step is to express the Heisenberg Hamiltonian in terms of ladder operators, ˆa, ˆa+ where [ˆa, ˆa+] = 1. Note that these are not the same operators used in the phonon derivation. Using the Holstein-Primakoff transformation the annihilation and creation operators become
Sˆi+= ˆSiz+ i ˆSiy =√ 2S
s
1 − ˆa+i ˆai
2S aˆi (18)
Sˆi−= ˆSix− i ˆSiy =√ 2Sˆa+i
s
1 − ˆa+i ˆai
2S (19)
when expressed in spherical coordinates[2]. Using the bosonic commutation relation and equation19the spin along the ˆz direction is
Sˆz= S − ˆa+ˆa (20)
This shows that the projection onto ˆz becomes the equlibrium minus a de- viation, see figure 1. Now the discrete Fouirer transform can be applied to ˆa and ˆa+ from real space to reciprocal space coordinates, assuming the spins are on a lattice structure with periodic boundary conditions, similar to the phonon derivation
ˆ ai=
r1 N
X
k
ˆbke−iksa (21)
a+i = r1
N X
k
ˆb+keiksa (22)
where k is the crystal momentum, N is the normalization constant and bk, b+k is the ladder operators in k space, or in other words the creation and annihilation operators for magnons. Using this transformation on equation18, 19and20and expanding the square roots in equation 18and19and assuming small deviations from equilibrium the Hamiltonian can be simplified to
H =ˆ 1
2J N zS2− 2µBB0S + ˆH0magnon+ O(1) (23) where O(1) is the order of accuracy that comes from the taylor expansion, z is the number of nearest neighbors and the Hamiltonian ˆH0magnonis the magnon spin in z direction,
Hˆ0magnon=X
k
ˆb+kˆbkωk~ (24)
where ωk = −2J zS(1 − γk) + 2µbHˆ0 is the dispersion relation, for magnons, γk =z1P
δeikδis the dispersion function and ˆb+kˆbk= ˆNk is the number operator representing the total number of magnons in mode k. Equation24can now be expressed similarly to equation14in its simplest form
Hˆ0magnon=X
k
Nˆk~ωk. (25)
Note that < ˆH >6=< ˆH0magnon>. The way of solving the spin interaction is similar to the phonon case. Using a transformation to linearizing the complex
problem and attaining magnons as a result. Even if the equations themselves are different the method still grants the same results. However in reality phonon and magnons interact with each other because the excange constant J dependes on the distance between the atoms.
3.4 The Coupling Of Phonon and Magnons
The phonon and magnons in a crystal structure depend on each other meaning that they are coupled. Applying a perturbation somewhere in the system and observing how the energy propagate throughout it, is dependent also on the crystal structure and material properties.
Figure 2: Illustrations of phonon-magnon interaction.
In figure 2the process is illustrated by blocks, phonons and balls, magnons connected with springs. When an energy is applied as in the figure2, the energy being in the form of motion, one can imagine how the motion distributes over the system. This is the fundamental idea of coupling between phonon and magnons.
3.5 The Boltzmann Equation
The Boltzmann equation is used in many areas of physics and in this project it will be used to describe the heat transport of phonon and magnons. The Boltzmann equation,
∂f
∂t = −dr dt
∂f
∂r −dp dt
∂f
∂p−∂f
∂t|collision+∂f
∂t|source, (26) describes an one particle density f moving and interacting with other particles or the surroundings. On the right side of equation26the first term describes the ballistic transport of the particle in time, the second part describes the change of momentum of the particle due to an external field, the third term represents the different interactions with other particles and the final term represent a source
term that generates new particles[4]. In case of an one dimensional problem, assuming the external force is 0, constant velocity and with no source term, the equation is reduced to
∂f
∂t + v∂f
∂x = ∂f
∂t|collision. (27)
In order to solve this equation the different collisions need to be taken into consideration. The first considered collision is when a phonon decays into a magnon and the second is when a magnon decays into a phonon. This is when only one phonon and one magnon is considered however in reality there are many more possible collisions that can occur. Considering a phonon mode (k, n) and magnon mode (q, p) where k, q are the wave vectors and n,p are polarization’s of respective quasi particle[4]. Each phonon and magnon mode have a specific group velocity and energy. For phonon: (vk,n, ωk,n) and for magnon: (cq,p, q,p).
The density of phonons and magnons described by the Boltzmann equation can be solved together
∂f (x, t)k,n
∂t + vk,n
∂f (x, t)k,n
∂x = ∂f (x, t)k,n
∂t |collision (28)
∂b(x, t)q,p
∂t + cq,p
∂b(x, t)q,p
∂x = ∂b(x, t)q,p
∂t |collision (29)
for phonon mode (k, n) and magnon mode (q, p). Phonon and magnon distributions at equilibrium are given by a Bose- Einstein distribution[4].
fk,nEq = 1 e
ωk,n
kB TEq − 1
bEqk,n= 1 e
k,n
kB TEq − 1
. (30)
Equation 30 depends on the temperature TEq which is space independent and kB which is the Boltzmann constant. When applying an external tempera- ture to an insulated one dimensional rod the overall temperature of the system will increase. The distribution will change with time before obtaining a new homogeneous temperature. It is essential to know the distribution of these par- ticles for every time step to solve the Boltzmann equation[5]. The collisions term is of high interest, describing the decay of quasi particles for a specific collision process. Applying it to a magnon - phonon coupling where a phonon can decay into a magnon and vice verse. All possible outcomes from processes of annihilation and creation of a phonon and magnon for the modes (k, n) and (q, p) can be described by,
∂fk,n
∂t |collision = γ1(fk,n+ 1)bq,p− γ2fk,n(bq,p+ 1) = 0 (31) where γ1 = τ1
k→q and γ2 = τ 1
q→k is inverse relaxation time for the different decays. The relaxation time τ is derived for one phonon and two magnon coupling, where the possible processes are a phonon colliding with a magnon to
create another magnon or two magnons decaying into a phonon and vice verse, by using the Fermi Golden rule[4].
1
τk→q =2π
~ D
k|H0|qE
2
δ(ωk− ωf). (32)
H0 is the Hamiltonian perturbation from equilibrium, k is the initial state while q is the final state. If the condition of energy conservation is not met the decay does not occur[6].
3.5.1 Equilibrium Temperature for Phonons and Magnons
Under equilibrium conditions equation 31 is balanced and therefore equal to zero. Inserting equation30into28and29the theoretical value for phonon and magnon equilibrium temperature can be calculated
kBTpTmln(γ1 γ2
) = Tpk,n− Tmωk,n (33) where Tp is the equilibrium phonon temperature and Tm is the equilibrium magnon temperature. To calculate the numerical value of the equilibrium tem- perature equation28and29can be rewritten to
Tp= ωk,n
kBln(h1
p+ 1) (34)
Tm= k,n kBln(h1
m + 1) (35)
where hp is the phonon distribution at equilibrium and hm is the magnon distribution at equilibrium. hp and hm are extracted from the simulation. If the calculated values from 34 and 35 satisfy equation 33 the numerical value matches the theoretical value.
3.6 Boundary and initial conditions for solving the Boltz- mann equation
In this project the Boltzmann equation will be solved along an one dimensional rod with length L, with reflective boundaries. This means that if a particle with velocity vk collides with a boundary the velocity becomes −vk. Initially, the system will be at thermal equlibrium, and then a thermal excitation δT will be applied. As a consequence phonons and magnons will be created, δfk,n(x, t) and δbk,n(x, t), respectively. The phonons and magnons will then propagate, decay and interact with each other. The computational approach to solve this will be explained in the Computational Method section.
4 Computational Method
The computational method used to simulate the phonon-magnon interaction is simple in theory. However it can become computationally heavy and difficult to program. This is the reason for using C + + which can handle big calculations effectively and make it possible to calculate the results in a reasonable amount of time. To analyze and plot the data Matlab was used because of the plot functions that Matlab has are simple and practical to use. The methods and approximations used to solve the system of coupled differential equations are described in the following sections.
4.1 Splitting Algorithm
In order to numerically solve equation27one can use the Strang operator split- ting algorithm[7]. Instead of looking at the transportation and collisions at the same time the equation can be separated into two different parts that happen after each other. First the particle is transported a time step ∆t2 and then the collision term is considerd over a time step ∆t[7]. The reason for different time step sizes for the transportation and collision term is because the collision term is more demanding. The Boltzmann equation can be written as
∂fk,n(x, t)
∂t = dfk,n(x, t)
dt |transport+dfk,n(x, t)
dt |collision (36) and utilizing the splitting algorithm the following system of equations are attained
∂fk,n(1)(x, t)
∂t =dfk,n(x, t)
dt |transport (37)
∂fk,n(x, t)
∂t =dfk,n(1)(x, t)
dt |collision (38)
where fk,n(1)(x, t) is the solution to the transport term[7]. The solution to 38 will be input in equation 37 and the process will repeat. For ∆t → 0 the splitting algorithm leads to the exact solution. Therefore it is relevant to work with small time steps.
To solve equation37and38the particle in cell approach will be used[7]. This method is used when simulating a continuous function over discrete intervals by creating pseudo particles. These pseudo particles are not real particles but helpful to simulate the behaviour of a continuous function in a computer[7]. By using a homogeneous distribution between the discrete points in the x plane and distributing pseudo phonons, or pseudo magnons, randomly across these homogeneous distribution it closely resembles a continuous function, given the pseudo particles are many. Every pseudo particle will have a position xk,n(t)
and a velocity vk,n. Reconstructing the function at a certain time is done with the following formula
fk,n(x, t) = C
N (t)
X
i=1
δ(x − xk’,n0(t))δ(k − k0) (39) where C is the normalization constant and N (t) is the number of pseudo particles at time t. Equation 39 basically counts how many pseudo particles there are between each discretization point and normalizes that quantity.
4.1.1 The Transport Term
Equation37describes a distribution function traveling with constant velocity.
This equation can be numerically solved with xk,n(t + dt) = xk,n(t) + vk,n
∆t
2 (40)
where xk,n(t) is a pseudo phonon at a certain position in time, vk,n is the group velocity for a certain mode k and ∆t2 is the time step size. It is also important to note that the Courant-Friedrich-Levy (CFL) condition
vmax∆t < ∆x (41)
for the maximum velocity vmax must be fulfilled[7]. This means that ∆t cannot be chosen independently of the discretization ∆x. It is practical to first choose a ∆x, all the velocities and then choose a ∆t so that the CFL condition is fulfilled.
4.1.2 The Collision Term
Equation38 describes the lifetime of a phonon or a magnon as a consiquence of the phonon-magnon scattering process. In the specific case studied in this project we consider the process of one phonon plus one magnon decaying into a magnon. This leads to different life times for the phonon and magnon modes.
Although the lifetimes for these processes have been computed, for sake of sim- plicity we will consider that the lifetimes are those of one phonon into one magnon or vice versa. Each process has a lifetime τk for every phonon mode and assuming the decay is exponential and dependent on the time step t the probability to decay becomes
P = 1 − e−∆tτk. (42)
If ∆t becomes large P → 1 which means that the probability for decay is almost 100%[7]. To simulate this decay a number between 0 and 1 will be generated using an uniform distribution. If the randomly generated number is smaller then P the decay occurs. If for instance a phonon pseudo particle decays it will create a magnon pseudo particle in the same place and vice versa.
If the randomly generated number is bigger than P the particle does not get removed. In case of one phonon one magnon the decay is only dependent on the random number however when more magnons or phonons are considerd the position needs to be taken into consideration.
This part of the code is very demanding because the collision term has to be calculated for every pseudo particle at every time step. Due to the many iterations required to simulate the transprot problem, the simulation time is very long.
4.1.3 Ratio of Minimum and Maximum Values
To investigate how fast the system reaches equilibrium the following equations will be used
αphonon=maxphonon(t) − minphonon(t)
maxphonon(t = 0) (43)
αmagnon=maxmagnon(t) − minmagnon(t)
maxmagnon(t = 0) (44)
where maxphonon(t), maxmagnon(t) is the maximum phonon and magnon population at time t and minphonon(t), minmagnon(t) is the minimum phonon and magnon population at time t. This ratio trends towards 0 when the system tend towards equilibrium.
5 Results
5.1 Numerical Values For Constants
In order to investigate the relaxation time of the system some numerical values are required. These values were calculated for NiO in the paper Microscopic Theory of Ultrafast Out-of-Equilibrium Dynamics in Magnetic Insulators. Un- raveling the Magnon-Phonon Coupling[6]. The in parameters that were kept constant during all simulations were the time step ∆t = 0.001ps, the length of the system Ltot = 4nm, the Boltzmann constant used for the distribution kb = 1.38064852 · 10−23J K−1,the discretisation in x, ∆x = 4 · 10−4, the number of pseudo particles Npseudo= 106the temperature peak of δT (x), Tmax= 300K, the width of δT (x), Twidth = 0.5. In this simulation δT (x) was shaped like a isosceles triangle with height Tmaxand width Twidth.
Table 1: Numerical values for relaxation time τk, group velocity vkand energies
k for phonon and magnons for a specific mode k. The values are calculated in the paper[6].
Decay time τk(ps) Group Velocity vk(m/s) Energy k(J )
Phonon 19.16 37.92 · 10−1 0.25 · 10−20
Magnon 5.46 200.000 · 10−1 1.56 · 10−20
5.2 Simulation With Initial Excited Phonon and Magnon
Utilizing the values from table 1 to excite both a phonon and a magnon and plotting the results for 1ps, 10ps, 70ps, 100ps, the stabilization time and equi- librium temperature can be approximated. Figures3, 4, 5, 6, 8 will show the phonon and magnon distributions at the respective times.
Figure 3: Phonon and magnon distribution using parameters from table1after 1ps. The initial phonon and magnon distribution was excited with δT (x) with peak at 300K and width of 0.5nm. The length of the system is given in nm
Figure 4: Phonon and magnon distribution using parameters from table1after 10ps. The initial phonon and magnon distribution was excited with δT (x) with peak at 300K and width of 0.5nm.The length of the system is given in nm
Figure 5: Phonon and magnon distribution using parameters from table1after 70ps. The initial phonon and magnon distribution was excited with δT (x) with peak at 300K and width of 0.5nm.The length of the system is given in nm
Figure 6: Phonon and magnon distribution using parameters from table1after 80ps. The initial phonon and magnon distribution was excited with δT (x) with peak at 300K and width of 0.5nm.The length of the system is given in nm
Figure 7: Phonon and magnon distribution using parameters from table1after 85ps. The initial phonon and magnon distribution was excited with δT (x) with peak at 300K and width of 0.5nm.The length of the system is given in nm
Figure 8: Phonon and magnon distribution using parameters from table1after 100ps. The initial phonon and magnon distribution was excited with δT (x) with peak at 300K and width of 0.5nm.The length of the system is given in nm
From figure 3, 4, 5, 6, 7 and 8 the maximum and minimum values of the phonon and magnon distributions were calculated. These values were calculated using Matlabs max and min function.
Table 2: Table of numerical values for maximum and minimum value for the phonon and magnon population at 1, 10, 70, 80, 85, and 100ps. The difference ratio of the maximum and minimum values over the initial maximum population is also displayed using equation43and44.
Time (ps) Magnon max value Phonon max value Magnon min value Phonon min value αmagnon αphonon
1 0.0018 7.87 ∗ 10−4 0 0 1 1
10 3.79 · 10−4 5.58 · 10−4 3 · 10−6 2 · 10−6 0.21 071
70 4.4 · 10−5 1.53 · 10−4 8 · 10−6 3.8 · 10−5 0.02 0.15
80 4.4 · 10−5 1.26 · 10−4 6 · 10−6 4.4 · 10−5 0.02 0.10
85 4.4 · 10−5 1.24 · 10−4 6 · 10−6 4 · 10−5 0.02 0.11
100 4 · 10−5 1.24 · 10−4 4 · 10−6 4.6 · 10−5 0.02 0.1
Using equation34and35the phonon temperature at equilibrium is 19.32K and the magnon temperature at equilibrium is 105.36K. The input for hp and hmwere calculadet using the average between minimum and maximum at 100ps.
The error of equation33was calculated with an error of order 10−19.
