Time-Domain Spectroscopy Tools

A small library for high-precision material parameter extraction from time-domain THz signals, written in Wolfram Language.

*Developed at Augsburg University, Germany*

__Github [page](https://github.com/JerryI/wl-tds-tools)__

Basics of THz spectroscopy

EM waves propagate through our sample and a reference (empty hole).

delta-pulseDetectorE(t)E(t)referencesample

Using gated detection, we sample electric field as a function of time. Using both signals we can reconstruct the responce function of the medium

$$

E(t) \\rightarrow E(\\omega), \\qquad \\hat{t} = \\frac{E(\\omega)^{sam}}{E(\\omega)^{ref}}, \\qquad Absorption \\approx - 2 ln(|\\hat{t}|)

$$

Installation

Using [LPM](https://github.com/JerryI/wl-localpackages)

PacletRepositories[{ Github -> "https://github.com/JerryI/wl-tds-tools" -> "master" }] Get["JerryI`TDSTools`Trace`"]; Get["JerryI`TDSTools`Transmission`"]; Get["JerryI`TDSTools`Material`"];

Let's get started

Firstly import our sample and reference data

(*BB[*)(*asumming you imported it somehow*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) ref = NotebookStore["ref"]; sam = NotebookStore["sam"]; ListLinePlot[{sam, ref}, PlotRange->Full, AxesLabel->{"t, ps", "I, nA"}] Create a time-trace object. We always enter used units explicitly

sam = TDTrace[QuantityArray[sam, {"Picoseconds", 1}]]; ref = TDTrace[QuantityArray[ref, {"Picoseconds", 1}]] List all properties

sam["Properties"] // TableForm Check the power spectrum and the range with a high signal-to-noise ratio

ListLinePlot[sam["PowerSpectrum"], PlotRange->Full] sam["FDCI"]

Reconstructing transmission function

On the next step we can derive an initial information about the transmission based on the sample and reference time-traces

tr = TransmissionObject[sam, ref, "Thickness"->Quantity[0.584994, "Millimeters"]] In the `TransmissionObject` decoration's preview, you can observe the actual transmission. A key challenge in digital signal processing is managing the phase of the signal. To address this, phase unwrapping is employed, which transforms a discontinuous phase signal into a continuous one.

There is a special function for it, but one have to play with a treshold

utr = TransmissionUnwrap[tr, "Basic", "PhaseThreshold"->5.3] Let us compare those two

ListLinePlot[{tr["Phase"], utr["Phase"]}, PlotLegends->{"Wrapped", "Unwrapped"}] However, you might need to tweak the settings, when we start to extract material parameters.

Material Parameters

This the main goal of all optical spectroscopy: use EM waves to probe the responce function of the material to extract all useful properties. It is sufficient to define $n$ and $\\kappa$ or their derivatives.

In the transmission spectroscopy we assume our crystal to be a slab taking into account all internal reflections. For the simplicity we can ignore them firstly

m1 = MaterialParameters[utr, "FabryPerotCancellation"->False] One can see already on the preview, that our refractive index is off. This rabid increate towards the low frequencies is a sign of a wrong phase offset. We can inspect on a regular plot

ListLinePlot[{ MaterialParameters[ Append[utr, "PhaseShift"->0] , "FabryPerotCancellation"->False]["n"] , MaterialParameters[ Append[utr, "PhaseShift"->-2] , "FabryPerotCancellation"->False]["n"] }, PlotRange->{{20,100}, {3.5,5.5}}, PlotLegends->{"PhaseShift = 0", "PhaseShift = -2"}] From here we one can guess that `-2` is a correct multipler. The next step will be to get rid of those wiggles we have from the reflections inside the medium

ListLinePlot[m1["k"], PlotRange->{{10,100}, {-0.03,0.2}}]

Fabry-Perot Cancellation

In general this is not a correct name (but quite recognizable). We have to deconvolve the material parameters from the signature of our sample shape. This can be done automatically using built-in iterative methods

{fp, nofp} = { MaterialParameters[ Append[utr, "PhaseShift"->-2] , "FabryPerotCancellation"->False] , MaterialParameters[ Append[utr, "PhaseShift"->-2] , "FabryPerotCancellation"->True] }; ListLinePlot[{fp["k"], nofp["k"]}, PlotRange->{{10,100}, {-0.03,0.2}}, PlotLegends->{"FP", "no FP"}] The success of this procedure strongly depends on phase, thickness accuracy and the correct ratio between sample and reference surface area.

Sometimes one have to vary the thickness, gain parameters to archive best results. Therefore it comes handy to have a criteria to judge the quality of the chosen combo

#["FPReduction"] &/@ {fp, nofp} This is done by making fourier transformation of the region of `"k"` with the highest signal-to-noise ratio, and then intergation over it. The wiggles usually give a sharp peak, which helps to estimate the magnitude of reduction effect.

Thickness and Gain correction

We can correct those parameters using `FPReduction` values as a criteria. The error mostly comes form the inaccuracy of thickness measurements, misaligment of reference/sample

map = Table[With[{ t = Append[utr, {"PhaseShift"->-2, "Gain"->gain, "Thickness"->Quantity[thickness,"Millimeters"]}] }, {thickness, gain, t} ], {thickness, 0.585-0.1, 0.585+0.1, 0.02}, {gain, 0.8, 1.2, 0.1}]; map = Flatten[map, 1]; map[[All,3]] = #["FPReduction"] &/@ MaterialParameters[map[[All,3]], "Target"->"CPU"]; ListContourPlot[map, PlotLegends->Automatic] Or if __GPU is available__, change the `Target` to `"GPU"`

map = Table[With[{ t = Append[utr, {"PhaseShift"->-2, "Gain"->gain, "Thickness"->Quantity[thickness,"Millimeters"]}] }, {thickness, gain, t} ], {thickness, 0.585-0.15, 0.585+0.15, 0.01}, {gain, 0.7, 1.2, 0.02}]; map = Flatten[map, 1]; map[[All,3]] = #["FPReduction"] &/@ MaterialParameters[map[[All,3]], "Target"->"GPU"]; ListContourPlot[map, PlotLegends->Automatic] An optimized pair can be found

optimized = MaximalBy[map, Last] // First Now let us plot the spectrum with corrected parameters

With[{m = MaterialParameters[Append[utr, { "Gain"->optimized[[2]], "PhaseShift"->-2, "Thickness" -> Quantity[optimized[[1]],"Millimeters"] }]]}, ListLinePlot[{m["Raw \[Alpha]"], m["\[Alpha]"]}, PlotRange->{{10,100}, {-5,100}}, PlotLegends->{"Raw", "Processed"}] ]

Live calculations

Interactive widget

Module[{ material, plotA, plotB }, {EventHandler[InputGroup[<| "Phase" -> InputRange[-2,2,1, "Label"->"Phase shift"], "Thickness" -> InputRange[0.5, 0.6, 0.005, "Label"->"Thickness"] |>], Function[result, material = MaterialParameters[ Append[utr, {"PhaseShift"->result["Phase"], "Thickness"->Quantity[result["Thickness"], "Millimeters"]}] , "FabryPerotCancellation"->True, "Target"->"GPU"]; {plotA, plotB} = QuantityMagnitude /@ {material["Best Raw \[Alpha]"], material["Best \[Alpha]"]}; ]] // EventFire, Graphics[{ColorData[97][1], Line[plotA//Offload], ColorData[97][2], Line[plotB//Offload]}, Frame->True, Axes->True] } // Row ]