Getting started

To get started we shall assemble a classical two-temperature model(TTM) in 0D for gold. For how to extend this to a 1D simulation and include thermal transport see Two Temperature Model. We shall use Unitful to make handling units easier. Unitful quantities can be supplied to any parameter in LightMatter.jl and the code will convert it to the correct unit system. For more on units see Unit Management

Creating our systems

The 3 main systems that comprise a Two-Temperature Model are an electronic temperature, a phononic temperature and a laser that drives the electronic system. The equation of motion for the TTM in 0D is given below:

\[ C_\text{el}\frac{\partial T_\text{el}(t)}{\partial t} = -g(T_\text{el} + T_\text{ph}) + S(t) \\ C_\text{ph}\frac{\partial T_\text{ph}(t)}{\partial t} = g(T_\text{el} + T_\text{ph})\]

Here, $S(t)$ is the equation for the laser, $C_\text{x}$ is the heatcapacity of the respective thermal bath and g is the electron-phonon coupling parameter.

In LightMatter.jl we need to create both thermal baths above as well as define an equation for our laser. This is all handed by 'build' functions. We shall build the most simple of each of the three systems in turn. First, the electronic thermal bath.

    Tel = build_ElectronicTemperature(Enabled=true, Electron_PhononCoupling=true, ElectronicHeatCapacity=:linear,
        ElectronPhononCouplingValue=:constant, γ = 62.9u"J/m^3/K^2", g = 0.3e17u"J/s/m^3/K")

ElectronicTemperature(true, false, true, false, :linear, :constant, 3.92590920783582e-7, 0.0, 0.0, 0.0, 1.872452722338229e-7)

This first system shows us the prototypical setup of building a system in LightMatter.jl. First, we turn the system on with 'Enabled = true'. Next we decide what approximations we want to use, in this case we are using a constant electron-phonon coupling value and a linear electronic heat capacity defined as $C_\text{el} = \gamma T_\text{el}$. Finally we then add the required values as Unitful quantities so we don't have to worry about units. For the phonon bath,

    Tph = build_PhononicTemperature(Enabled=true, Electron_PhononCoupling=true, PhononicHeatCapacity=:constant,
    Cph = 2.49e6u"J/m^3/K")

PhononicTemperature(true, false, true, false, :constant, 0.0, 0.0, 0.0155413575954073, 0.0)

The same principle applies here as for the electronic bath, but now we don't need to provide details about the coupling as the phonon bath assumes that there is an electronic system to couple to and the parameters are located there. We decided on a constant value of the phonon heat capacity and this is the only value we need to provide. Last but not least is the laser.

    Las = build_Laser(envelope=:Gaussian, FWHM=50.0u"fs", ϕ=10.0u"J/m^2", Transport=:optical, ϵ=16.3e-9u"m")

Laser(:Gaussian, :optical, 50.0, 62.41509074460764, 0.0, 0.0, 16.3, 0.0, 0.0)

Here we have chosen a gaussian envelope for the laser profile. To see the other options see Lasers. We have provided a full-width half-maximum (FWHM) and a fluence (ϕ). The last two settings describe how the laser interacts with the material. :optical states that the interaction depends on the absorption coefficient (1/ϵ). ϵ is then the inverse of the absorption coefficient a.k.a the penetration depth.

Once we define all of our systems we combine them into the Simulation.

    Sim = build_Simulation(electronictemperature = Tel, phononictemperature = Tph, laser=Las)

LightMatter.Simulation(ElectronicTemperature(true, false, true, false, :linear, :constant, 3.92590920783582e-7, 0.0, 0.0, 0.0, 1.872452722338229e-7), PhononicTemperature(true, false, true, false, :constant, 0.0, 0.0, 0.0155413575954073, 0.0), AthermalElectrons(false, false, false, false, false, :FLT, :constant, :unity, :constant, 0.0, 0.0, 0.0, [NaN]), ElectronicDistribution(false, false, 1.0, 5.685621837291226), PhononicDistribution(false, false, 0.0, 3-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Flat()) with element type Float64:
 Ratios.SimpleRatio{Int64}(4, 1)
 Ratios.SimpleRatio{Int64}(5, 1)
 Ratios.SimpleRatio{Int64}(6, 1), 0.0), Structure(false, false, 1, 3-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Int64}, Gridded(Linear())), Flat()) with element type Float64:
 Ratios.SimpleRatio{Int64}(4, 1)
 Ratios.SimpleRatio{Int64}(5, 1)
 Ratios.SimpleRatio{Int64}(6, 1), DataInterpolations.AkimaInterpolation{Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64}[DataInterpolations.AkimaInterpolation{Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64}([1, 2, 3], [4, 5, 6], Float64[], [1.0, 1.0, 1.0], [0.0, 0.0], [0.0, 0.0], DataInterpolations.ExtrapolationType.Constant, DataInterpolations.ExtrapolationType.Constant, FindFirstFunctions.Guesser{Vector{Int64}}([4, 5, 6], Base.RefValue{Int64}(1), true), false, false), DataInterpolations.AkimaInterpolation{Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64}([1, 2, 3], [4, 5, 6], Float64[], [1.0, 1.0, 1.0], [0.0, 0.0], [0.0, 0.0], DataInterpolations.ExtrapolationType.Constant, DataInterpolations.ExtrapolationType.Constant, FindFirstFunctions.Guesser{Vector{Int64}}([4, 5, 6], Base.RefValue{Int64}(1), true), false, false)], [-10.0, -9.99, -9.98, -9.97, -9.96, -9.95, -9.94, -9.93, -9.92, -9.91  …  9.91, 9.92, 9.93, 9.94, 9.95, 9.96, 9.97, 9.98, 9.99, 10.0], LightMatter.Dimension(1, [0.0], 1.0, 0.0), LightMatter.TotalFields(LightMatter.Fields([0.0, 0.0, 0.0], [0.0, 0.0, 0.0]), LightMatter.Fields([0.0, 0.0, 0.0], [0.0, 0.0, 0.0]))), Laser(:Gaussian, :optical, 50.0, 62.41509074460764, 0.0, 0.0, 16.3, 0.0, 0.0), DensityMatrix(false, Matrix{Complex}[[0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im], [0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im], [0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im]], LightMatter.Fields(Expr[:((.+)((sqrt.(((2 * sim.laser.ϕ * sqrt(4 * log(2))) .* 1.0) ./ (Constants.c * Constants.ϵ0 * sim.laser.FWHM * sim.laser.n * sqrt(pi))) .* exp((-2 * log(2) * t ^ 2) / sim.laser.FWHM ^ 2)) .* cos.(((sim.laser.hv * 1) / (Constants.ħ * 2 * pi)) * t), 0.0)), :((.+)(0.0 + 0.0, 0.0)), :((.+)(0.0 + 0.0, 0.0))], Expr[:((.+)(0.0 + 0.0, 0.0)), :((.+)(Expr[:((sqrt.(((2 * sim.laser.ϕ * sqrt(4 * log(2))) .* 1.0) ./ (Constants.c * Constants.ϵ0 * sim.laser.FWHM * sim.laser.n * sqrt(pi))) .* exp((-2 * log(2) * t ^ 2) / sim.laser.FWHM ^ 2)) .* cos.(((sim.laser.hv * 1) / (Constants.ħ * 2 * pi)) * t)), :(0.0 + 0.0), :(0.0 + 0.0)] ./ Constants.c, 0.0)), :((.+)(0.0 + 0.0, 0.0))]), Complex[0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im]))

Using the Simulation

Now that we have the simulation we can go straight to running it if we desire. But before that we can have a look at the equations of motion that LightMatter.jl has constructed for our simulation. To do this we run,

    eom = LightMatter.function_builder(Sim)

Dict{String, Union{Expr, Vector{Expr}}} with 2 entries:
  "Tph" => :((0.0 + -1 * (-(sim.electronictemperature.g) * (Tel - Tph))) ./ sim…
  "Tel" => :((((sqrt((4 * log(2)) / pi) / sim.laser.FWHM) * exp((-4 * log(2) * …

This returns us a dictionary for each of the equations of motion. In this case the keys are "Tel" and "Tph". We can have a look at the equation of motion for the electronic thermal bath.

    eom["Tel"]

:((((sqrt((4 * log(2)) / pi) / sim.laser.FWHM) * exp((-4 * log(2) * t ^ 2) / sim.laser.FWHM ^ 2)) * (1 ./ sim.laser.ϵ) * ((1 .- sim.laser.R) .* sim.laser.ϕ) + -(sim.electronictemperature.g) * (Tel - Tph)) / (sim.electronictemperature.γ * Tel))

This equation looks like what we would expect. We can see the beginning is the equation for a normalised gaussian laser. This is then followed by the subtraction of the electron-phonon coupling parameter multiplied by the temperature difference of the thermal baths and this is all divided by the expression for the heat capacity that we expected.

We can now run a simulation by defining a couple initial temperatures and a time span,

    initialtemps=Dict("Tel"=>300.0,"Tph"=>300.0)
    tspan=(-150.0, 1000.0)

(-150.0, 1000.0)

Finally, we can now run the simulation. We use the OrdinaryDiffEq package for performing the time integration and so all arguments used within a standard solve call can be used here. You can also use the alg keyword to define any time integration algorithm supported by OrdinaryDiffEq.

    sol = run_simulation(Sim, initialtemps, tspan, saveat=1.0, abstol=1e-10, reltol=1e-10)

retcode: Success
Interpolation: 1st order linear
t: 1151-element Vector{Float64}:
 -150.0
 -149.0
 -148.0
 -147.0
 -146.0
 -145.0
 -144.0
 -143.0
 -142.0
 -141.0
    ⋮
  992.0
  993.0
  994.0
  995.0
  996.0
  997.0
  998.0
  999.0
 1000.0
u: 1151-element Vector{RecursiveArrayTools.NamedArrayPartition{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, @NamedTuple{Tph::Int64, Tel::Int64}}}:
 (Tph = [300.0], Tel = [300.0])
 (Tph = [300.00000000000006], Tel = [300.000000010539])
 (Tph = [300.0000000000002], Tel = [300.00000002518357])
 (Tph = [300.0000000000007], Tel = [300.0000000454816])
 (Tph = [300.00000000000136], Tel = [300.00000007361575])
 (Tph = [300.0000000000025], Tel = [300.0000001126737])
 (Tph = [300.00000000000415], Tel = [300.0000001664157])
 (Tph = [300.0000000000066], Tel = [300.0000002398222])
 (Tph = [300.00000000001006], Tel = [300.0000003406232])
 (Tph = [300.00000000001495], Tel = [300.0000004793324])
 ⋮
 (Tph = [346.93990297645223], Tel = [3985.0558443600876])
 (Tph = [346.9837328115791], Tel = [3984.620424881217])
 (Tph = [347.02755687271133], Tel = [3984.1850151867657])
 (Tph = [347.07137516003667], Tel = [3983.7496152780773])
 (Tph = [347.1151876737423], Tel = [3983.3142251564987])
 (Tph = [347.1589944140158], Tel = [3982.8788448233754])
 (Tph = [347.20279538104467], Tel = [3982.4434742800527])
 (Tph = [347.24659057501634], Tel = [3982.008113527879])
 (Tph = [347.29037999611853], Tel = [3981.572762568199])

Due to using RecursiveArrayTools the solution output can be quite tricky to use. To make this easier for instant plotting / processing you can call the following to output a dictionary of the solution of each system as well as the save times as Array's.

    results = LightMatter.seperate_results(sol, initialtemps, Sim)

Dict{String, Union{Float64, AbstractArray}} with 4 entries:
  "Tph"   => [300.0; 300.0; … ; 347.247; 347.29;;]
  "noe"   => [40.0;;]
  "Tel"   => [300.0; 300.0; … ; 3982.01; 3981.57;;]
  "times" => [-150.0, -149.0, -148.0, -147.0, -146.0, -145.0, -144.0, -143.0, -…

The last thing to complete our first simulation is to save the results to a file. This can be handled by post_production. This will save the results and Simulation struct to a HDF5 file as well as can perform some automated post processing. To see what can be done check out Outputting.

    post_production(sol, "Test.hdf5", initialtemps, [:ThermalFermiDistribution], Sim)

What's next?

Now that we've covered the basics of performing a TTM, we're ready to explore all the capabilities of LightMatter.jl All the systems follow the same patterns of enabling, defining approximations and adding parameters but with a couple of extra notes when you start to perform more complicated simulations. Feel free to check out the Tutorials or Systems sections to found out more. The tutorials have the structure as follows: Theory of the method implemented, Approximations that you can use, Performing the method in LightMatter.jl.