import numpy as np
import escape as esc
esc.require('0.9.1')
from escape.utils import show
import escape.utils.mdb as mdb
x
# Gamma emissions under grazing incidence and resonance-enhanced neutron standing waves in a thin film
X-ray or neutron standing waves can be generated above a mirror surface under conditions of total external reflection by incident and reflected beams.
These beams have equal amplitude and superpose to form standing waves inside the sample. With the proper sample morphology, i.e. layer thickness close to standing waves period, the amplitude of the standing waves can be
many-times amplified, causing emission of X-ray fluorescense in the case of X-rays or gammas in the case of neutrons by a layer
of "heavier" material or stronger absorber. In this notebook we demonstrate how `ESCAPE` can be used for the modelling of Gamma emissions together with neutrons specular reflectivity.
As a basis of our model and obtained results we use the following publication:
`Grazing incidence prompt gamma emissions and resonance-enhanced neutron standing waves in a thin film`
Huai Zhang, P. D. Gallagher, S. K. Satija, R. M. Lindstrom, R. L. Paul, T. P. Russell, P. Lambooy, and E. J. Kramer
Phys. Rev. Lett. 72, 3044
https://doi.org/10.1103/PhysRevLett.72.3044
## Theoretical background
The calculation of the intensity of Gamma emissions starts with the the Schrödinger equation which solution is the wave function in each layer given by
$\psi_m(z) = T_m exp(ik^{(m)}_z(z-z_m))+R_m exp(-ik^{(m)}_z(z-z_m))$,
where $z_m$ corresponds to the the upper interface of the $m$-th layer, $T_m$ and $R_m$ are the trasmission and reflection Fresnel coefficients.
Continuity of the wave functions and their derivatives at each interface allows to find Fresnel coefficients for each layer using so-called `matrix formalism`.
The analysis of conservation of probability in the case of absorption, i.e. non-zero imaginary part of the scattering potential leads to the relation for the neutron capture rate
of a single isotope element given by
$\tau=c\frac{4\pi}{k_z}\int\left|\psi(z)\right|^2\rho_{im}dz$,
where $k_z$ is the transverse wave vector component, $\psi(z)$ - neutron wave function, $\rho_{im}$ - imaginary part of scattering length density,
$c$ - normalization constant. The inegration is performed over the layer $Z$ - boundaries.
Further details can be found in the original publication or in the third edition of `Quantum Mechahics` by Leonard I. Schiff.
X-ray or neutron standing waves can be generated above a mirror surface under conditions of total external reflection by incident and reflected beams.
These beams have equal amplitude and superpose to form standing waves inside the sample. With the proper sample morphology, i.e. layer thickness close to standing waves period, the amplitude of the standing waves can be
many-times amplified, causing emission of X-ray fluorescense in the case of X-rays or gammas in the case of neutrons by a layer
of "heavier" material or stronger absorber. In this notebook we demonstrate how ESCAPE
can be used for the modelling of Gamma emissions together with neutrons specular reflectivity.
As a basis of our model and obtained results we use the following publication:
Grazing incidence prompt gamma emissions and resonance-enhanced neutron standing waves in a thin film
Huai Zhang, P. D. Gallagher, S. K. Satija, R. M. Lindstrom, R. L. Paul, T. P. Russell, P. Lambooy, and E. J. Kramer Phys. Rev. Lett. 72, 3044
https://doi.org/10.1103/PhysRevLett.72.3044
The calculation of the intensity of Gamma emissions starts with the the Schrödinger equation which solution is the wave function in each layer given by
,
where corresponds to the the upper interface of the -th layer, and are the trasmission and reflection Fresnel coefficients.
Continuity of the wave functions and their derivatives at each interface allows to find Fresnel coefficients for each layer using so-called matrix formalism
.
The analysis of conservation of probability in the case of absorption, i.e. non-zero imaginary part of the scattering potential leads to the relation for the neutron capture rate of a single isotope element given by
,
where is the transverse wave vector component, - neutron wave function, - imaginary part of scattering length density, - normalization constant. The inegration is performed over the layer - boundaries.
Further details can be found in the original publication or in the third edition of Quantum Mechahics
by Leonard I. Schiff.
xxxxxxxxxx
## Sample description
According to the publication the sample consists of the following layers starting from Si substrate: Si/Ni(600 nm)/PVCH(72nm)/Gd(5nm)/PVCH(72nm)/dPS(25nm).
First we should create materials instances. For the polymer layers we use the mass density as fit parameter. The imaginary part of the scattering legth density is very small for these layers
and it's influence on the reflectivity is negligible compared to Gd layer. Plus, we reduce the number of fit parameters.
### Materials
According to the publication the sample consists of the following layers starting from Si substrate: Si/Ni(600 nm)/PVCH(72nm)/Gd(5nm)/PVCH(72nm)/dPS(25nm).
First we should create materials instances. For the polymer layers we use the mass density as fit parameter. The imaginary part of the scattering legth density is very small for these layers and it's influence on the reflectivity is negligible compared to Gd layer. Plus, we reduce the number of fit parameters.
dPS = mdb.amorphous("C8D8",
density=esc.par(name="C8D8", value=0.95, units="g/cm^3"),
source=mdb.Source.NEUTRONS) # deuterated polystyrene
PVCH = mdb.amorphous("C8H14",
density=esc.par(name="PVCH", value=0.9, units="g/cm^3"),
source=mdb.Source.NEUTRONS) # polyvinylcyclohexane
Gd = mdb.amorphous("Gd", source=mdb.Source.NEUTRONS)
Ni = mdb.amorphous("Ni", source=mdb.Source.NEUTRONS)
xxxxxxxxxx
### Layers
Initial thickness and roughness values for each layer are equal to the values provided in the publication or estimated from the provided SLD profile. In our model roughness from interfaces is included as in the standard matrix method statistically averaged intensity over the interfaces roughness. In the publication interface roughness was modelled by smearing the abrupt change in the scattering potential at each interface with a hyperbolic tangent function. It make sense if a very large interface roughness is expected, when a standard statistical averaging of thickness changes perpendicular to the interface plane
Initial thickness and roughness values for each layer are equal to the values provided in the publication or estimated from the provided SLD profile. In our model roughness from interfaces is included as in the standard matrix method statistically averaged intensity over the interfaces roughness. In the publication interface roughness was modelled by smearing the abrupt change in the scattering potential at each interface with a hyperbolic tangent function. It make sense if a very large interface roughness is expected, when a standard statistical averaging of thickness changes perpendicular to the interface plane
angstr=u"\u212B"
dps_layer = esc.layer("Layer: dPS", dPS,
thkn=esc.par(value=240, units=angstr), rough=esc.par(value=15, units=angstr))
pvch1_layer = esc.layer("Layer: PVCH1", PVCH,
thkn=esc.par(value=737, units=angstr),
rough=esc.par(value=14, userlim=[0, 20], units=angstr))
pvch2_layer = esc.layer("Layer: PVCH2", PVCH,
thkn=esc.par(value=741, units=angstr),
rough=esc.par(value=14, userlim=[0, 20], units=angstr))
gd_layer = esc.layer("Layer: Gd", Gd,
thkn=esc.par(value=84, userlim=[40, 120], units=angstr),
rough=esc.par(value=14, userlim=[0, 20], units=angstr))
substr = esc.substrate("Substrate: Ni", Ni,
rough=esc.par(value=30, userlim=[0, 100], units=angstr))
sample = esc.multilayer("dPS/PVCH/Gd/PVCH/Ni/Si", esc.air("Air"), substr)
sample.add(dps_layer)
sample.add(pvch1_layer)
sample.add(gd_layer)
sample.add(pvch2_layer)
Z=esc.var("Z")
prof=esc.profile_sld0("Profile SLD0", Z, sample)
show(esc.real(prof), coordinates=np.linspace(-100, 2500, 5000), xlabel=u"Z[\u212B]", ylabel=u"SLDRe[\u212B^-2]", xlog=False)
xxxxxxxxxx
## Specular reflectivity
Gamma-ray intensities and neutron reflectivity data were measured simulataneously keeping the illuminated area of the sample constant at $~35\times40\;mm^2$ with the beam resolution $\Delta k_z/k_z\approx0.058$. The latter is used for averaging of specular reflectivity over the resolution. The FWHM of the resolution function is a function of $k_z$ in this case.
Gamma-ray intensities and neutron reflectivity data were measured simulataneously keeping the illuminated area of the sample constant at with the beam resolution . The latter is used for averaging of specular reflectivity over the resolution. The FWHM of the resolution function is a function of in this case.
Kz0 = esc.var("Kz0") #variable after integration over resolution function
fwhm = Kz0*0.058 # FWHM is a functor
Kz = esc.var("Kz") # variable of ideal reflectivity functor (no average)
Rmod = esc.specrefl("Refl", Kz*2, sample, formalism="matrix")
I0 = esc.par("I0", 0.7)
B = esc.par("Bckgr", 1, scale=1e-5)
#averaging
Rmod_av = I0*esc.average_normal( Rmod, fwhm, Kz, Kz0, maxiter=300, numstd=10, epsabs=1e-8)+B
#after averaging R is a function of Kz0
#print(Rmod_av.variables)
#The experimental data has been obtained with the data scanner from the original publication.
# As a result the error information has been lost.
#Below we added small errors to every point according to the Poisson distribution.
#The value of incident beam intensity (1e7) is used to get realistic errorbar values.
Kz_data, R_data = np.loadtxt("data/reflectivity.csv", unpack=True)
errors=np.sqrt(R_data*1e7)/np.max(R_data*1e7)
data=esc.data("Reflectivity", Kz_data, R_data, errors, copy=True)
mdl=esc.model("Reflectivity Model", Rmod_av, data, "data", "q4")
xxxxxxxxxx
## Gamma emissions
In `ESCAPE` we introduced `reftrans` package which implements matrix formalism and allows to calculate Fresnel coefficients inside particular layer of a multilayer sample as well as the wave function. In the next cell we demonstrate how to use them for the calculation of Gamma emissions. The neutron capture rate should be finally multiplied by a constant C, representing the counting efficiency. According to the publication, C was one of the fit parameters together with structural parameters.
In ESCAPE
we introduced reftrans
package which implements matrix formalism and allows to calculate Fresnel coefficients inside particular layer of a multilayer sample as well as the wave function. In the next cell we demonstrate how to use them for the calculation of Gamma emissions. The neutron capture rate should be finally multiplied by a constant C, representing the counting efficiency. According to the publication, C was one of the fit parameters together with structural parameters.
#Gd layer has index 3
gd_thkn=gd_layer.parameter(2)
gd_sldi=gd_layer.parameter(1)
Rm = esc.reftrans_rm("Rm", Kz, sample, 3) # reflection coefficient inside Gd layer
Tm = esc.reftrans_tm("Tm", Kz, sample, 3) # transmission coefficient inside Gd layer
Pm = esc.reftrans_kz("Pm", Kz, sample, 3) # transversal component of the wave-vector inside Gd layer
xxxxxxxxxx
The imaginary part of the scattering potential $u_{im}^{(m)}$ is constant inside the layer and the solution of the integral $\int_{Z_m}\left|\psi_m\right|^2dz$ is straightforward:
The imaginary part of the scattering potential is constant inside the layer and the solution of the integral is straightforward:
a = esc.real(Pm)
b = esc.imag(Pm)
cRm = esc.conj(Rm)
cTm = esc.conj(Tm)
tm=gd_thkn
psi_int_2 = (-esc.norm(Tm)/(2*b)*(1-esc.exp(2*b*tm))
+esc.norm(Rm)/(2*b)*(1-esc.exp(-2*b*tm))
+Tm*cRm/(2*1j*a)*(1-esc.exp(-2*1j*a*tm))
-cTm*Rm/(2*1j*a)*(1-esc.exp(+2*1j*a*tm))
)
x
# imaginary part of the scattering potential is given by gd_sldi*4*pi
g_rate = esc.real(psi_int_2)*esc.abs(gd_sldi)*4*np.pi/Kz
Ig0 = esc.par("Ig0", 1.3e-2)
#averaging over the resolution function
g_rate_av = Ig0*esc.average_normal( g_rate, fwhm, Kz, Kz0, maxiter=50, numstd=3)
Kzg_data, G_data = np.loadtxt("data/gamma.csv", unpack=True)
errors=np.zeros(shape=G_data.shape)
datag=esc.data("Gamma", Kzg_data, G_data, errors, copy=True)
mdlg=esc.model("Gamma Model", g_rate_av, datag, "none", "none")
#we create now optimizer and perform fitting with the user interface
opt = esc.diffevol("", [mdl, mdlg], maxiter=25)
show(opt, ylog=[True, False], xlabel=["Kz [1/\u212B]", "Kz [1/\u212B]"], ylabel=["Reflectivity", "Gamma rate"])