THz Time-domain modeling

Using a single Lorentz-oscillator model, we model the response of a material with variable thickness.

Download Notebook

Get the complete code

When electromagnetic radiation propagates through the material, the transmitted radiation field can be described as:

E(\omega)^{\prime} = t(\omega) E(\omega)^{(source)}

where $t$ is a complex transmission function describing the material response. For a relatively thick material, neglecting Fresnel losses, you can relate it to material parameters such as the refractive index $n$ :

t \approx e^{- i \omega n / (c L)}

where $L$ is the plate thickness. For non-magnetic systems, the material response can be described by $\varepsilon \approx n^2$ , the complex dielectric function. A common model for the dielectric response is the Lorentz oscillator, which describes the motion of bound charges, ionic motion, or electronic excitations:

\varepsilon(\omega) = \varepsilon_{\infty} + \frac{\omega_p^2}{\omega_0^2 - \omega^2 + i \gamma \omega}

which has resonance frequency $\omega_0$ , strength $\omega_p^2$ , and damping $\gamma$ (losses). This expression is an approximation valid over a limited frequency range; contributions from higher-lying excitations are approximated by the constant background $\varepsilon_{\infty}$ , which can be set to 1.

Let's define the dielectric function f for a single Lorentz oscillator:

ClearAll[f,n,e,c,tds,tds2line];

f[ω_, eps_, ω0_, ωp_, γ_] := eps + (*FB[*)(((*SpB[*)Power[ωp(*|*),(*|*)2](*]SpB*))(*,*)/(*,*)((*SpB[*)Power[ω0(*|*),(*|*)2](*]SpB*) - (*SpB[*)Power[ω(*|*),(*|*)2](*]SpB*) + I γ ω))(*]FB*);

n[ω_, params__] := (*SqB[*)Sqrt[f[ω, params]](*]SqB*)

The real (phase) and imaginary (loss) parts of n look like this:

Plot[n[x, 9.8, 28.66, 15, 1.5] // ReIm // Evaluate, {x,28.66 - 20, 28.66 + 20}, 
  FrameLabel->{"wavenumber (cm^{-1})"}, Frame->True
]

(*VB[*)(FrontEndRef["8d1a6108-bdf6-4de9-85d2-dbddb2e17679"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKW6QYJpoZGljoJqWkmemapKRa6lqYphjppiSlpCQZpRqam5lbAgCMCRYQ"*)(*]VB*)

This gives the material response in the frequency domain. We are interested in the time domain, which is what is measured in experiments.

Now substitute it into the transmission function e, passing all arguments:


e[ω_, L_, params__] := Exp[- (*FB[*)((I 2π (*SpB[*)Power[10(*|*),(*|*)12](*]SpB*) ω n[ω, params] L)(*,*)/(*,*)(c cm2THz))(*]FB*) ] /.{c -> 3.0 (*SpB[*)Power[10(*|*),(*|*)10](*]SpB*), cm2THz -> 33.4}

In THz spectroscopy, a short broadband pulse is typically used to probe the material; the response is measured as the oscillating electric field versus time. Here we model that pulse as a delta function in the frequency domain:

\delta(\omega) = 1

To switch to the time domain we use Fourier together with some signal processing to obtain the real time-domain signal:

tds[params__] := (Join[#, Reverse[#]] &@ Table[e[ω, params], {ω, 0,100,0.1}]) // Fourier
tds2line[params__] := {Range[Length[#]], #} &@ (Drop[#, -Length[#]/2//Ceiling] &@ tds[params] // Re) // Transpose

Finally, plot the resulting radiation (projection of the E field) as a function of time, keeping the material thickness as a variable parameter:

Manipulate[Graphics[{
  ColorData[97][1],
    Line[tds2line[thickness/10.0, 9.8, 28., 15., 0.05]], Black,
    Text[StringTemplate["d = `` mm"][thickness], {300, 2.5}]
  },
  
  PlotRange->{{0,400}, {-3,3}}, 
  Axes->True, Frame->True, FrameLabel->{"time, ps", "Electric field, nA"}
], {{thickness, 0.5, "d (mm)"}, 0.01, 1, 0.02},
  ContinuousAction -> True,
  Appearance->None
]

(*VB[*)(FrontEndRef["a2e7c659-bc03-46ce-b662-6ed193b3ae55"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxqlmiebmVrqJiUbGOuamCWn6iaZmRnpmqWmGFoaJxknppqaAgCJSxXU"*)(*]VB*)

Even with fixed material parameters, the response changes with thickness due to propagation effects, producing nonlinear beatings. This phenomenon - polariton beatings or interference - occurs when the light field hybridizes with the Lorentz-oscillator mode.

Applications

Low-lying phonons, magnons, excitons, crystal-field transitions and all kind of species, which have a pole in Green's function. Here is a fitted time-spectrum of $Co_{2} Mo_{3} O_{8}$ antiferromagnet with a strong magnon mode oscillating at $\sim 40~cm^{-1}$ :

(*VB[*)(Legended[ToExpression[FrontEndRef["7496de5b-d37f-4368-873a-d8b4ad7cbb9d"], InputForm], Placed[LineLegend[{Directive[PointSize[0.006944444444444445], RGBColor[0.24, 0.6, 0.8], AbsoluteThickness[2], Opacity[0.9]], Directive[PointSize[0.006944444444444445], RGBColor[0.95, 0.627, 0.1425], AbsoluteThickness[2], Opacity[0.9]]}, {"Exp.", "Model"}, LegendMarkers -> {{False, Automatic}, {False, Automatic}}, Joined -> {True, True}, LabelStyle -> {}, LegendLayout -> {"Column", {Automatic, 5}}], {0.6, 0.3}, Identity]])(*,*)(*"1:eJylUssuREEQHe9HSLBGIrEdCzPMWExuPIaQmXhcse++XZeO1j26+2LE1kewtfMF1jY2WLDwAWyESPyBfoTJndiIWpxUV1dXn6pTY1hsxG2ZTEb1G9iicLgAkZBICxl2mUgFtoGTybjVpvQaKBNq7mxi3GJjQwYWpeC6zEn5CKJEI8wgHDfhQn5mmsAUzpJcIc7mc9PFbLGQQ1lSxHlEChHGM8QXbjewkZhn3dYBRFY5q7vopkzA8+s0sMZQBCTu/CZToRw8w0adClU6to7qMbBAJUSaHoBna0NrgnId0mOQIzf7wzf7J4H/wP29NDcvmJDyavT0bf3qPpA5Z8+BPD+z9hr4QoMGZrESLNGwuUOjXQ5KUUvC39vZrdZQRHVd3t1aew/+zyp29hHIy8+HKh54CWSp9+miVrr+D6vU5ELrlI9qE2GHcaqCAEtL5PbEz7yK5C5I1TT61EnZKouIKXBNziZa7CFNo79muZ+t7CtmTNCk9s+iNDYm/dCtCsLAQl1nEGd+5+tS+366q6C6SJraCS0FI0Wyx5sopHhT2086obFH3nkMnLLLBLg2anwBHW/rpQ=="*)(*]VB*)

However there is some background presented and a few weaker modes at higher frequnecies, therfore this fit is not really "perfect", but is good enough for our purposes.

Since we work in the linear optics regime, nothing can stop you from adding more terms to f function.

Bonus

Source code of the preview image:


{Gray, {Directive["Roughness"->0.0], Cuboid[{Pi, -10.0/3, -10/2}, {1.4Pi,10.0/2,10.0/3}]}, Directive["Emissive"->Red, "EmissiveIntensity"->10],Red, Cases[Plot[5Exp[-(*SpB[*)Power[(x)(*|*),(*|*)2](*]SpB*)]Cos[10 x], {x,-2Pi,Pi}, PlotRange->Full][[1]], _Line, Infinity] /. {Line[arr_] :> Tube[Map[Function[xy, {xy[[1]], 0.3, xy[[2]]}], arr], 0.1]}};
Graphics3D[%,  "Renderer"->"PathTracing"]