Ammonia Maser

Introduction: The (Subscript[NH, 3])molecule and its potential

The ammonia molecule ((Subscript[NH, 3])) consists of a pyramid with the N at the top or at the bottom:

Off[General::spell1]
pts         =  {{0, 0, 0}, {1.2, 0, 0}, {0.5, (Sqrt[5/4]), 
            0}, {0.5, (Sqrt[5/4])/2, 2/(Sqrt[8])}};
siz         =  { .09, 0.09, 
          0.09,  .15};
col         =  {1, 1, 1,  .6};
points  = 
        Table[{Hue[col[[i]]], PointSize[siz[[i]]], 
            Point[pts[[i]]]}, {i, 1, 4}];
lines     =  
        Map[Line, 
          Flatten[Table[
              Table[{pts[[i]], pts[[j]]}, {j, i + 1, 4}], {i, 1, 3}], 
            1]];
dl             =  {Dashing[{0.05, 0.05}], 
          Thickness[ .01], 
          Line[{{pts[[4, 1]], pts[[4, 2]], -pts[[4, 3]]}, 
              pts[[4]]}]};

    
pp = 
        Graphics3D[{PointSize[0.08], points, lines, dl}, 
          PlotRange -> {{- .1, 1.2}, {0, 1.5}, {- .6,  .8}}, 
          Boxed -> False, Axes -> False, 
          ViewPoint -> {1.786,  -2.497,  
              1.422}];
Show[pp];

The three H-atoms can be above or below the nitrogen. Both configurations represent a minimum in the potential energy and the H-atoms can swing from one position to the other. This transition is used in an ammonia maser to produce microwaves. Alternatively, they can be taken as fixed and the nitrogen is then moving in a potential with two minima separated by a finite potential barrier.

Save plot to ps-file

SetDirectory["~/mathem"];
Export["ammonia.eps", pp, "EPS"];

Potential

The simplified potential which is going to be treated, is shown here:

Definition of figure

$TextStyle = {FontFamily -> "Times", FontSize -> 16, 
          FontWeight -> "Bold"};
g1 = 
        Graphics[{Thickness[ .02], 
            Line[{{-1, 100}, {-1, 0}, {- .25, 0}, {- .25, 
                  1}, { .25, 1}, { .25, 0}, {1, 0}, {1, 100}}]}, 
          PlotRange -> {{-1.5, 1.5}, {- .5, 2}}];
\(g2 = \ 
      Graphics[{
			Text[\*""(Subscript[V, 0]) →"", {\(-1.3\), 
              1}], 
			Text["0 →", {\(-1.3\), 
              0.05}], 
			Text["b", { .25, \(- .2\)}], 
			Text["-b", {\(- .25\), \(- .2\)}], 
			Text["-a", {\(-1\), \(- .2\)}], 
			Text["a", {1, \(- .2\)}], 
			Text["←  Δ  →", {\(- .6\),  .2}], 
			Text["←  Δ  →", { .6,  .2}], 
			Text["↑ ∞", {1.2, 1.9}]
				}, 
        PlotRange -> {{\(-1.5\), 1.5}, {\(- .5\), 2}}]; 
    Show[{g1, g2}];\)

Meaning of Parameters:

Δ  =  a - b
β  =  (b/Δ)

Save plot to ps-file

SetDirectory["~/mathem"];
Export["nh3_pot.eps", pp, "EPS"];

Parameters

\(Clear["Global`*"]; \  (*Clear\ all\ previous\ global\ settings\ from\ currently\ open\ sessions*) 
Needs["Graphics`Shapes`"];\)

Parameters of the potential (see introduction):

a0 = 1.3; b0 = 0.3; V00 = 18;

Number of timesteps for the simulation:

Nt0 = 50;

Schroedinger equation

The complete potential for the N-atom has several parts:

V=∞		for	 	|x| >a
V=0		for		|x| < b
V = (Subscript[V, 0])	   for		b <|x|< a

This gives the time independent Schroedinger equation as:

ψ(x)=0						for			|x| ≥a
-((ℏ^2)/(2 m)) ((d^2)/(dx^2)) Ψ(x)  +  (Subscript[V, 0])x Ψ(x)  =  E Ψ(x)		for			|x| < b
-((ℏ^2)/(2 m)) ((d^2)/(dx^2)) Ψ(x)  =  E Ψ(x)	  		             	for			b <|x|< a

Introducing k =  ((Sqrt[2 m E])/ℏ) and (Subscript[k, V]) =  ((Sqrt[2 m (Subscript[V, 0])])/ℏ), these equations can be expressed as:

ψ(x)=0						for			|x| >a
((d^2)/(dx^2)) Ψ(x)    =  - ((k^2) - (Subscript[k, V]^2)) Ψ(x)			for			|x| < b
 ((d^2)/(dx^2)) Ψ(x)  =  -(k^2) Ψ(x)	  		            	for			b <|x|< a

Here, we are only interested in solutions with (k^2)<(Subscript[k, V]^2).

Parts of eigenfunctions

Solution for -a < x ≤ -b, using  ψ(-a)=0:

(Subscript[L, 1])[k_] = 
      Block[{a, b}, 
        Simplify[
            DSolve[{(Subscript[ψ, 1])''[x] + (k^2) (Subscript[ψ, 1])[x] == 
                  0, (Subscript[ψ, 1])[-a] == 0}, (Subscript[ψ, 1])[x], x], 
            k ∈ Reals] /. {C[2] -> A , 
            C[1] -> B}]

Solution for b ≤ x ≤ a, using  ψ(a)=0:

(Subscript[L, 3])[k_] = 
      Simplify[DSolve[{(Subscript[ψ, 3])''[x] + (k^2) (Subscript[ψ, 3])[x] == 
                0, (Subscript[ψ, 3])[a] == 0}, (Subscript[ψ, 3])[x], x], 
          k ∈ Reals] /. {C[2] -> C , 
          C[1] -> D}

Solution for -b ≤ x ≤ b:

(Subscript[L, 2])[k_] = 
      Simplify[DSolve[{(Subscript[ψ, 2])''[
                    x] + ((k^2) - (kv^2)) (Subscript[ψ, 2])[x] == 
                0}, (Subscript[ψ, 2])[x], x], {kv ∈ Reals, 
            k ∈ Reals}] /. {C[2] -> F, 
          C[1] -> G}

Eigenvalues and Eigenfunctions

As the potential is symmetric with respect to the y-axis, the eigenfunctions have to be either symmetric or antisymmetric. This can be exploited to simplify the coefficients and gives:

Symmetric:		C = - A	,	G = F
Antisymmetric:     	 C = A		  ,	  G = - F

\(\(ψ\_1[k_, x_] = \(ψ\_1[x] /. \(L\_1[k]\)[\([1]\)]\) /. 
          A -> A\&~/Sec[a\ k];\)
 (*D[ψ\_1[k, x], 
        x]*) \)
\(\(ψ\_\(3  S\)[k_, x_] = \(ψ\_3[x] /. \(L\_3[k]\)[\([1]\)]\) /. 
          C -> \(-A\&~\)/Sec[a\ k];\)\)
\(\(ψ\_\(3  A\)[k_, x_] = \(ψ\_3[x] /. \(L\_3[k]\)[\([1]\)]\) /. 
          C -> A\&~/Sec[a\ k];\)\)
\(\(ψ\_\(2  S\)[k_, x_] = 
        FullSimplify[\(\(ψ\_2[x] /. \(L\_2[k]\)[\([1]\)]\) /. 
              G -> F\) /. 
            F -> F\&~/2];\)
 (*D[ψ\_\(2  A\)[k, 
          x], x]*) \)
\(\(ψ\_\(2  A\)[k_, x_] = 
        FullSimplify[\(\(ψ\_2[x] /. \(L\_2[k]\)[\([1]\)]\) /. 
              G -> \(-F\)\) /. 
            F -> F\&~/2];\)\)
\(\(Print[
        TableForm[{{\*""(Subscript[ψ, 1])="", ψ\_1[k, 
                x], " ", \*""(Subscript[ψ, 1])'="", 
              D[ψ\_1[k, x], 
                x]}, {\*""(Subscript[ψ, 2  S])="", ψ\_\(2  S\)[
                k, x], " ", \*""(Subscript[ψ, 2  S])'="", 
              D[ψ\_\(2  S\)[k, x], 
                x]}, \ 
{\*""(Subscript[ψ, 2  A])="", ψ\_\(2  A\)[k, 
                x], " ", \*""(Subscript[ψ, 2  A])'="", 
              D[ψ\_\(2  A\)[k, x], 
                x]}, \ 
{\*""(Subscript[ψ, 3  S])="", ψ\_\(3  S\)[k, 
                x], " ", \*""(Subscript[ψ, 3  S])'="", 
              D[ψ\_\(3  S\)[k, x], 
                x]}, \ 
{\*""(Subscript[ψ, 3  A])="", ψ\_\(3  A\)[k, 
                x], " ", \*""(Subscript[ψ, 3  A])'="", 
              D[ψ\_\(3  A\)[k, x], x]}}]];\)\ \)

The wavefunction must be continuously differentiable, this condition yields:

(Subscript[ψ, 1])(-b)=(Subscript[ψ, 2  S])(-b) 	,	(Subscript[(ψ'), 1]) (-b)=(Subscript[(ψ'), 2  S]) (-b)		or:   	(Subscript[ψ, 1])(-b)=(Subscript[ψ, 2  A])(-b) 		,	(Subscript[(ψ'), 1]) (-b)=(Subscript[(ψ'), 2  A]) (-b)

Each pair of equations gives a system with the two independent variables  \(A\&~\)and \(F\&~\):

Symmetric:	      				             Sin [k (a - b)]\(A\&~\)  - \(Cosh[\(-b\)\ \@\(\(-k\^2\) + k\_V\%2\)]\ F\&~\) = 0
					\(k\ Cos[k\ \((a - b)\)]\ A\&~\) - \(Sinh[\(-b\) \@\(\(-k\^2\) + k\_V\%2\)]\ \(\@\(\(-k\^2\) + k\_V\%2\)\) 
        F\&~\)=0

Antisymmetric:			                                   Sin [k (a - b)]\(A\&~\) +Sinh[-b (Sqrt[-(k^2) + (Subscript[k, V]^2)])]\(F\&~\) = 0		
					\(k\ Cos[k\ \((a - b)\)]\ A\&~\) + \(Cosh[\(-b\) \@\(\(-k\^2\) + k\_V\%2\)]\ \(\@\(\(-k\^2\) + k\_V\%2\)\) 
        F\&~\)=0

These equations lead to the introduction of new, dimensionless quantities: kd = k (a-b), β = b / (a-b) and kvd = (Subscript[k, V]) (a-b). In order to give a non-trivial solution, the determinant of the above systems have to be 0. This gives us:

Symmetric:		tan(kd)=- \(TraditionalForm\`\_kd\/\@\(kvd\^2 - kd\^2\)\)coth[ β(Sqrt[(kvd^2) - (kd^2)])]
Antisymmetric:		tan(kd)=- \(TraditionalForm\`\_kd\/\@\(kvd\^2 - kd\^2\)\)tanh[ β(Sqrt[(kvd^2) - (kd^2)])]

kvd := (Sqrt[2  V0])(a - b);β := b/(a - b);

Eigenvalues

We solve these equations for | kd | < | kvd | approximately:

\(Block[{a, b, V0, nn, ϵ, k, k0, k1, n, 
        km}, 
a = a0; b = b0; 
      V0 = V00; 
ϵ\ \ \  = \  .00000000001; 
nn\  = \ Floor[\(2\ kvd\  + \ π\)\/\(2\ π\)]; 
fs[
          kd_]\  := \ Tan\ [kd] + 
          kd\/\@\(kvd\^2 - kd\^2\)\ Coth[\@\(kvd\^2 - kd\^2\)\ β]; 

      fa[kd_]\  := 
        Tan\ [kd] + 
          kd\/\@\(kvd\^2 - kd\^2\)\ Tanh[\@\(kvd\^2 - kd\^2\)\ β]; 
Plot[{Tan\ [
            kd], \(-\(kd\/\@\(kvd\^2 - kd\^2\)\)\)\ Coth[\@\(kvd\^2 - kd\^2\)\ β], \(-\(kd\/\@\(kvd\^2 - kd\^2\)\)\)\ Tanh[\@\(kvd\^2 - kd\^2\)\ β]}, {kd, 0, kvd}, PlotRange -> {{0, kvd}, {\(-5\), 5}}, 
        AxesLabel -> {"kd", ""}, 
        PlotStyle -> {Hue[ .6], Hue[1], Hue[0.3]}, 
        TextStyle -> {FontWeight -> 
              Plain}]; 
ks\  = 
        Flatten[\ 
          Table[\ k1\  = \ Min[n\ π, kvd] - ϵ; \ 
            k0\  = \ n\ π - π/2 + ϵ; \ 
            km\  = \ \((k0 + k1)\)/2; 
			If[fs[k0]*fs[k1] < 0, 
              k /. FindRoot[
                  fs[k] == 0, {k, km, k0, 
                    k1}], {}]
\ \ \ \ \ \ \ \ \ \ , {n, 1, 
              nn}]]; 
ka\  = \ 
        Flatten[Table[\ k1\  = \ Min[n\ π, kvd] - ϵ; \ 
            k0\  = \ n\ π - π/2 + ϵ; \ 
            km\  = \ \((k0 + k1)\)/2; 
			If[fa[k0]*fa[k1] < 0, 
              k /. FindRoot[
                  fa[k] == 0, {k, km, k0, 
                    k1}], {}]
, {n, 1, 
              nn}]];
]\)

The intersections of the blue Tangent-curves with the red and green curves give the possible values of kd for the symmetric and antisymmetric wavefunctions.

Plot of eigenfunctions

For the obtained k, we can choose \(TraditionalForm\`A\&~\)and \(TraditionalForm\`F\&~\) to define continuos  wave-functions , the constants are finally determined by the normalization.

\(\(φs[i_, x_] := 
        If[x < \(-1\) - β, 0, 
          If[x < \ \(-β\), \ 
            As\&~[\([i]\)]\ \ Sin[ks[\([i]\)] \((1 + β + x)\)], 
            If[x\  < \ β, 
              Fs\&~[\([i]\)]\ Cosh[\@\(\(kvd\^2\)\(-\)\(\(\(ks[\([i]\)]\)\(\ \)\)\^2\)\(\ \)\)\ x], 
              If[x < 1 + β, 
                As\&~[\([i]\)]\ Sin[ks[\([i]\)] \((1 + β - x)\)], 
                0]]]];\)\)
\(\(φa[i_, x_] := 
        If[x < \(-1\) - β, 0, 
          If[x < \ \(-β\), \ 
            Aa\&~[\([i]\)]\ \ Sin[ka[\([i]\)] \((1 + β + x)\)], 
            If[x\  < \ β, 
              Fa\&~[\([i]\)]\ Sinh[\@\(\(kvd\^2\)\(-\)\(\(\(ka[\([i]\)]\)\(\ \)\)\^2\)\(\ \)\)\ x], 
              If[x < 1 + β, \(-Aa\&~[\([i]\)]\)\ Sin[
                    ka[\([i]\)] \((1 + β - x)\)], 0]]]];\)\)

\(\(Block[{a, b, V0, i, x}, 
b = b0; a = a0; 
        V0 = V00; 
As\&~ = 
          Cosh[\@\(\(kvd\^2\)\(-\)\(ks\^2\)\(\ \)\)\ β]; 
        Fs\&~ = \ Sin[ks]; 
Zs = 
          Table[NIntegrate[Abs[φs[i, x]]\^2, {x, \(-a\), a}], {i, 
              1, Length[ks]}]; 
As\&~ = 
          Cosh[\@\(\(kvd\^2\)\(-\)\(ks\^2\)\(\ \)\)\ β]/\@Zs; 
        Fs\&~ = \ Sin[ks]/\@Zs; 
Aa\&~ = 
          Sinh[\@\(\(kvd\^2\)\(-\)\(ka\^2\)\(\ \)\)\ β]; 
        Fa\&~ = \ \(-Sin[ka]\); 
Za = 
          Table[NIntegrate[Abs[φa[i, x]]\^2, {x, \(-a\), a}], {i, 
              1, Length[ka]}]; 
Aa\&~ = 
          Sinh[\@\(\(kvd\^2\)\(-\)\(ka\^2\)\(\ \)\)\ β]/\@Za; 
        Fa\&~ = \ \(-Sin[ka]\)/\@Za; 
psymm = 
          Plot[{φs[1, x], φa[1, x]}, {x, \(-a\), a}, 
            AxesLabel -> {\*""\!\(x\/\(a - b\)\)"", ""}, 
            PlotLabel -> "symmetric", 
            PlotStyle -> {Hue[1], Hue[0.65]}, 
            TextStyle -> {FontWeight -> 
                  Plain}]; 
pasymm = 
          Plot[{φs[2, x], φa[2, x]}, {x, \(-a\), a}, 
            AxesLabel -> {\*""\!\(x\/\(a - b\)\)"", ""}, 
            PlotLabel -> "antisymmetric", 
            PlotStyle -> {Hue[1], Hue[0.65]}, 
            TextStyle -> {FontWeight -> 
                  Plain}];
];\)\)

Store plots to ps-files

SetDirectory["~/mathem"];
Export["symmetric.eps", 
        psymm, "EPS"];
Export["antisymm.eps", pasymm, "EPS"];

Dynamics

Definition

\(\(\(
\)\(dyn[ca_, cs_, tmax_, Nt_] := 
        Block[{x, na = Length[ca], ns = Length[cs], j, p1, it, 
            dt = tmax/Nt}, 
F[x_, t_] := 
            Sum[Exp[I\ t\ \((ks[\([j]\)])\)\^2/
                        2]\ cs[\([j]\)]*φs[j, x], {j, 1, ns}] + 
              Sum[Exp[I\ t\ \((ka[\([j]\)])\)\^2/2] 
                  ca[\([j]\)]*φa[j, x], {j, 1, 
                  na}]; 
p1 = 
            Plot[Abs[F[x, 0]]\^2, {x, \(-a\), a}, 
              PlotRange -> {{\(-a\), a}, {0, 2}}, 
              DisplayFunction -> Identity]; 
Table[
            Show[{p1, 
Plot[
                  Abs[F[x, dt\ it]]\^2, {x, \(-a\), a}, 
                  PlotStyle -> {Hue[1]}, 
                  PlotRange -> {{\(-a\), a}, {0, 2}}, 
                  DisplayFunction -> Identity]}, 
              AxesLabel -> {"x", "ρ(x)"}, 
              TextStyle -> {FontWeight -> Plain}, 
              DisplayFunction -> $DisplayFunction], {it, 1, 
              Nt}]; 
	
          SelectionMove[EvaluationNotebook[], After, 
            EvaluationCell]; 
	
          SelectionMove[EvaluationNotebook[], Next, 
            CellGroup]; 
	
          FrontEndTokenExecute["OpenCloseGroup"];
];\)\)\)

This simulation shows the dynamics of a superposition of the first two eigenstates (i.e. of one symmetric and one antisymmetric eigenstate):

Block[{a = a0, b = b0, V0 = V00, ca, 
          cs}, 
  
ca = {1/N[(Sqrt[2])], 
            0}; 
cs = {1/N[(Sqrt[2])], 0}; 
dyn[
          ca, cs, 4 π/((((ka[[1]]))^2) - (((ks[[1]]))^2)), 
          Nt0]];

The period of this inter-well oscillation is (Subscript[τ, 1])=\(TraditionalForm\`\(4  π\)\/\(k\_\(1, S\)\%2 - k\_\(1, A\)\%2\)\):

\(\(\(
\)\(Print[\*""\!\(τ\_1\)= "", 
        4\ π/\((\((ka[\([1]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\)];\)\)\)

This simulation shows the dynamics of a superposition of the first two symmetric (i. e. the first and the third) eigenstates:

Block[{a = a0, b = b0, V0 = V00, ca, 
          cs}, 
ca = {0, 
            0}; 
cs = {1/N[(Sqrt[2])], 
            1/N[(Sqrt[2])]}; 
dyn[ca, cs, 
          4 π/((((ks[[2]]))^2) - (((ks[[1]]))^2)), 
          Nt0]];

The period of this intra-well oscillation is (Subscript[τ, 2])=((4  π)/((Subscript[k, 2, S]^2) - 
                      (Subscript[k, 1, S]^2)))≪(Subscript[τ, 1]).

Nt2 = 
        Round[2 ((((ks[[2]]))^2) - (((ks[[1]]))^2))/((((ka[[1]]))^2) - (((ks[[1]]))^2))];
\(\(Print[\*""(Subscript[τ, 2])= "", 
        4\ π/\((\((ks[\([2]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\), \*""       (Subscript[τ, 1])/(Subscript[τ, 2])= "", \((\((ks[\([2]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\)/\((\((ka[\([1]\)])\)\^2 - \((ks[\([1]\)])\)\^2)\)];\)\)

The timescale of these two motions can be seen in this simulation, which shows a superposition of the first three eigenstates:

● Caution: The number of timesteps is roughly 2(Subscript[τ, 1])/(Subscript[τ, 2]) which can be rather large!

Block[{a = a0, b = b0, V0 = V00, ca, 
          cs}, 
cs = {0.5, 
            0.5}; 
ca = {1/(Sqrt[2]), 0}; 
dyn[
          ca, cs, 4 π/((((ka[[1]]))^2) - (((ks[[1]]))^2)), 
          Nt2]];