In [17]:
import numpy as np
import escape as esc
esc.require('0.9.7')
from escape.utils.widgets import show

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.

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¶

In [18]:
dPS = esc.amorphous("DPS", formula="C8H[2]8",  density="0.95g/cm^3")
PVCH = esc.amorphous("PVCH", formula="C8H14",  density="0.9g/cm^3")
Gd = esc.amorphous("Gd", density=None, dwsld0re="1+-1", dwsld0im="1+-1")
Ni = esc.amorphous("Ni", density=None, dwsld0re="1+-1", dwsld0im="1+-1")
print(PVCH, dPS)

Name: PVCH                	Parameters number:          1
Parameter           	Value          	+-	Error     	Units     	Fixed
PVCH Density        	            0.9	+-	0         	g/cm^3    	    0
 
Name: DPS                 	Parameters number:          1
Parameter           	Value          	+-	Error     	Units     	Fixed
DPS Density         	           0.95	+-	0         	g/cm^3    	    0

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

In [19]:
dps_layer = esc.layer("Layer: dPS", dPS, thkn=esc.par(value=240, units=esc.angstr), rough=esc.par(value=15, units=esc.angstr))
pvch1_layer = esc.layer("Layer: PVCH1", PVCH, 
                        thkn=esc.par(value=737, units=esc.angstr),
                       rough=esc.par(value=14, userlim=[0, 20], units=esc.angstr))
pvch2_layer = esc.layer("Layer: PVCH2", PVCH,
                        thkn=esc.par(value=741, units=esc.angstr),
                        rough=esc.par(value=14, userlim=[0, 20], units=esc.angstr))
gd_layer = esc.layer("Layer: Gd", Gd,
                     thkn=esc.par(value=84, userlim=[40, 120], units=esc.angstr),
                     rough=esc.par(value=14, userlim=[0, 20], units=esc.angstr))
substr = esc.substrate("Substrate: Ni", Ni,
                       rough=esc.par(value=30, userlim=[0, 100], units=esc.angstr))

sample = esc.multilayer("dPS/PVCH/Gd/PVCH/Ni/Si", frgr=esc.air("Air"), bkgr=substr)

sample.add(dps_layer)
sample.add(pvch1_layer)
sample.add(gd_layer)
sample.add(pvch2_layer)

src = esc.neutrons(4.1, units="angstrom")
show(sample, source=src)

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.

In [20]:
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", source=src)


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

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 [21]:
#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, source=src) # reflection coefficient inside Gd layer
Tm = esc.reftrans_tm("Tm", Kz, sample, 3, source=src) # transmission coefficient inside Gd layer
Pm = esc.reftrans_kz("Pm", Kz, sample, 3, source=src) # transversal component of the wave-vector inside Gd layer

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:

In [23]:
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 [24]:
# 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", 5, scale=1e-8)

#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", "q4") 

show(mdlg, xlabel="Kz [1/\u212B]", ylabel="Gamma rate [n/sec]", title="Gamma emissions")

Wave-function functor¶

Instead of finding the analytical solution the user could use wave-function functor (reftrans_wf method) and numerical integration.

In [25]:
Z=esc.var("Z")
#Wave function is a functor of two variables Kz and Z

dps_thkn=dps_layer.parameter(1)
pvch1_thkn=pvch1_layer.parameter(1)

psi = esc.reftrans_wf("Psi", Kz, Z, sample, source=src)

#numerical integration of |psi|^2 inside the layer
ll = dps_thkn + pvch1_thkn # upper interface of Gd layer
ul = ll + gd_thkn          # lower interface

psi_int = esc.integral(esc.norm(psi), Z, ll, ul, maxiter=50)

g_rate = Ig0*psi_int*esc.abs(gd_sldi)*4*np.pi/Kz

#averaging over resolution function
g_rate_av = esc.average_normal( g_rate, fwhm, Kz, Kz0, maxiter=50, numstd=3)

mdlg=esc.model("Gamma Model", g_rate_av, datag, "none", "none") 

#Finally we get the same result, but calculation is slower
show(mdlg, xlabel="Kz [1/\u212B]", ylabel="Gamma rate [n/sec]", title="Gamma emissions")
In [ ]:
In [26]:
#we create now optimizer and perform fitting with the user interface
opt = esc.diffevol("", [mdl, mdlg], maxiter=25)
#opt()#remove this line to skip preliminary fit
show(opt, ylog=[True, False], xlabel=["Kz [1/\u212B]", "Kz [1/\u212B]"], ylabel=["Reflectivity", "Gamma rate"])

We can also build a 2D image representing probability density along the sample normal for the full range of Z and Kz arguments. Using this representation one can effectively search for the conditions of resonance inside layers of interest.

In [27]:
kz=np.linspace(0, 0.08, 100)
za=np.linspace(-800, 2600, 1500)

# create meshgrid of both coordinates
xv, yv=np.meshgrid(kz, za)
# create coordinates array of type [x0, y0, x1, y1,...]
coords=np.column_stack([xv.flatten(), yv.flatten()]).flatten()

show(esc.norm(psi), coordinates=coords, xlabel="Kz [1/\u212B]", ylabel="Z[\u212B]", plot_type="map", rows=1500, 
     coord_index=[0, 1], cblog=False, cbmin=1e-2)

The resonance Kz-positions of the resonance enhanced neutron standing waves for several modes can be seen in this figure. Dragging thickness and SLD sliders of Gd and capping layers it is possible to visualize a strong dependence of resonance conditions on structural parameters.

Using reduce method which allows to reduce number of variables substituing constants or parameters instead of variables, we can make slicing of the 2D image and have a look at the wave-function amplitude along resonance modes, for example, 4th and 5th as in the publication.

In [28]:
psi_n4=esc.reduce("",  esc.norm(psi), Kz, 0.00707)
psi_n5=esc.reduce("",  esc.norm(psi), Kz, 0.0089)
show(psi_n4, coordinates=np.linspace(-800, 2600, 5000), title="n=4 mode")
show(psi_n5, coordinates=np.linspace(-800, 2600, 5000), title="n=5 mode")

For the $n=4$, the position probability density given by $|\psi|^2$ is resonantly amplified in the PVCH film and for the $n=5$ the standing wave is suppressed by the absorption in the Gd layer. Due to the imperfect data and fit results you probably will not see this effect after fitting. The wave function is very sensitive to the density of the capping layer and thicknesses. Dragging the corresponding sliders the effect can be significantly increased keeping the model parameters in their physical ranges.

Roughness model limits¶

In this notebook we have demonstrated how to build the model of Gamma emissions together with specular neutron reflectivity. Qualitatively, our results are very close to the published results in the reference (Zhang et al, 1994). Quantitatively, there are expected discrepancies in the curves due to the lack of information and different model for the roughness.

In our calculations roughness has a statistical nature and the resulting reflectivity and gamma emission values are statistically averaged. In the matrix formalism the statistically averaged amplitudes cannot be found analytically and are approximated. This approximation is valid for the relatively small roughness values.

The probability plot for the whole sample can be used to check the validity of the used model. Wave-function as well as its derivative should be continuous at the interfaces. Negligible differences are acceptable, or you should probably use another approach, for example, gradient layers.

In [29]:
show(psi_n4, coordinates=np.linspace(1700, 1900, 100), title="n=4 mode")
In [ ]: