In [20]:

`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

${\psi}_{m}(z)={T}_{m}exp(i{k}_{z}^{(m)}(z-{z}_{m}))+{R}_{m}exp(-i{k}_{z}^{(m)}(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 {|\psi (z)|}^{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.

`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.

In [22]:

`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

In [23]:

`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 $\text{}35\times 40\phantom{\rule{thickmathspace}{0ex}}m{m}^{2}$ with the beam resolution $\mathrm{\Delta}{k}_{z}/{k}_{z}\approx 0.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.

In [24]:

`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") `

``

In [ ]:

``

`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.

In [25]:

`#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 ${u}_{im}^{(m)}$ is constant inside the layer and the solution of the integral ${\int}_{{Z}_{m}}{\left|{\psi}_{m}\right|}^{2}dz$ is straightforward:

In [26]:

`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))`

` )`

``

In [29]:

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") `

``

In [30]:

`#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"])`