Computing the S-matrix for an odd pair of delta functions.

Let’s compute the S-matrix for the odd delta potential,
        V(x) = g((ℏ^2)/(2m L))[δ(x+L)-δ(x-L)]
where g is a dimensionless number.  Note that the boundary condition at the deltas is thus, eg at x=+L,
        (∂_x)ψ((SuperPlus[L])) =(∂_x)ψ((SuperMinus[L]))-(g/L)ψ(L) 
Since the potential is constant away from the deltas, we can expand the solution in each region as the sum of plane waves,
        ψ = (({{A (e^( i k x)) +B (e^( -i k x))                       x<-L}, {({{F (e^( i k x))+G (e^( -i k x))        -L<x<+L}, {C (e^( i k x))+D (e^( -i k x))        +L<x             }})}})) 
Imposing the boundary conditions at each delta will give us four equations relating the six coefficients.  The S-matrix is the matrix that gives the outgoing coefficients when acting on the ingoing coefficients, ie
        (({{B}, {C}}))=S (({{A}, {D}}))
We'll deduce this by plugging our boundary conditions into mmatica, solving them with Solve, and extracting the coefficients.  

Before we get started, let' s declare all the variables we' ll be using to be real and positive -- when working with complex numbers in Mathematica, this greatly simplifies your life.

$Assumptions={β>0,k>0,L>0,g>0,η>0};

Thus prepared, let's input the equations representing the boundary conditions at the walls of the wells (ie, continuity of the wavefunction and its derivative).   Note that i’ve named the variables AA, BB, CC, DD, FF and GG, to avoid confusion with global variables, and defined a replacement rule to introduce the dimensionless momentum, η=2k L.

ψ1 = AA (E^(I k x))+BB (E^(-I k x));
ψ2 = FF (E^(I k x))+GG (E^(-I k x));
ψ3 = CC (E^(I k x))+DD (E^(-I k x));

EQs={ ψ1==ψ2/.{x->- L},      (∂_x)ψ2==(∂_x)ψ1+(g/L)ψ1/.{x->- L},
        ψ2==ψ3/.{x->L},         (∂_x)ψ3==(∂_x)ψ2-(g/L)ψ2/.{x->L}}//Simplify
ktoE={k->(η/(2L))};

To solve the equations, we use the command Solve.  The basic syntax of Solve is Solve[{eqn1, eqn2, eqn3...}, {var1, var2, ...}], where eqn1 is the first equation (be sure to use double equal signs == in your equations), var1 is the first variable you want to solve for, and so on.  We have four equations and six variables, so the best we can do is determine four in terms of the remaining two.  To get the S-matrix, we want the independent variables to be the ingoing components, AA and DD, so:

Sol = Collect[Solve[EQs/.ktoE,{FF,GG,BB,CC}][[1]]//Expand//FullSimplify,{AA,DD}];

Note that we’ve already replaced k with η.   We can now extract the S-matrix by pulling off the coefficients of AA nad DD in the solution for BB and CC.  Mathematica provides a nice way to do this in the form of the CoefficientArrays function.  Calls to this function take the form {b,S1,S2...}=CoefficientArrays[{eqn1,eqn2,...},{var1,var2,...}] where {var1,var2,...}  are the independent variables, b is the vector of constants in the equations, S1 is the matrix of terms linear in the variables, S2 is the matrix of terms quadratic in the variables, etc.  for us, the equations are linear, so we only need two terms:

{b,S}=CoefficientArrays[{BB+(BB/.Sol)==0,CC+(CC/.Sol)==0},{AA,DD}]//Normal//FullSimplify;
MatrixForm[S]

As a quick check, let’s verify that our S-matrix is Unitary:

MatrixForm[S\[ConjugateTranspose].S//ComplexExpand//Simplify]

Nice.

Now that we are confident of our S-matrix, let’s extract some info from the system.

First, note that while (Subscript[S, 12])= (Subscript[S, 21]), we find that (Subscript[S, 11])≠ (Subscript[S, 22]), which reflects that the potential does not respect Parity.  (Side question.  As you saw on the problem set, a potential which is even under Parity leads to an S-matrix which has (Subscript[S, 12])= (Subscript[S, 21]) and  (Subscript[S, 11])= (Subscript[S, 22]).  In the present case, with a potential that is odd under Parity, we find that (Subscript[S, 12])= (Subscript[S, 21]) and  (Subscript[S, 11])=(Subscript[S, 22])*.  Is this a general rule?  What does this tell you about how S, as a matrix, transforms under Parity?  What is the physical meaning of this difference?)

Next let’s look at the transmission probability and the scattering phase.

t=S[[2,1]];
T = (((Abs[t]//ComplexExpand))^2);  (* Do a little brute-force simplification... *)
T = (T//TrigExpand)/.{((Cos[η])^2)->1-((Sin[η])^2)}//Simplify
ϕ =ArcTan[( -((Im[t])/(Re[t]))//ComplexExpand)]//Simplify

Plotting the transmission probability and scattering phase make the results particularly obvious.  First the probability:

Manipulate[Plot[T/.{η->ηη,g->gg}, {ηη, 0, 5π}, PlotRange->{-0.05, 1.05}], {gg, 0.001, 10}]

Showing the expected resonances at η=N π.

Now the phase shift:

Manipulate[Plot[ϕ/.{η->ηη,g->gg}, {ηη, 0, 5π}, PlotRange->{-(π/2)-0.05,(π/2)+0.05}], {gg, 0.001, 10}]

Question: what can you read off this plot?

Now let’s lok for bound states, ie states with E<0.  As discussed in lecture, bound states correspond to poles of S: when we take k->i α, so that (e^( i k x))->(e^( -α x)), normalizability requires A,D=0 but B,C≠0, which is only possible if S has a pole.  So: if we take k->i α, what must be true such that S has a pole?  To get units straight, recall that η=2k L -> 2i α L = i β, so that β=2α L.

MatrixForm[Sb= S/.{η->I β}//Simplify]

Poles of S are zeroes of the denominator for which the numerator is non-zero,

Denominator[Sb[[1,2]]]==0//FullSimplify

That’s transcendental, so we may as well solve it graphically.  A particularly convenient way to do so is to write the equation as,
        (g^2)==((β^2)/(1-(e^(-2 β))))
and plot LHS and RHS together as a function of g.  Manipulate makes this easy:

Manipulate[Plot[{(g^2),((β^2)/(1-(E^(-2β))))},{β,0,5π},PlotRange->All],{g,0.001,10}]

We thus have a single bound state for all values of g.  Note:  It may be tempting to think that there’s a threshold bound state (β=0)  at g=0;  however, the numerators also vanish at  β=0, so this is not a pole of S, and correspondingly there is no bound state.  Instead, this is the usual zero-energy state of the free particle.

Note that, when we first looked for the bound states, we could also have taken k->-i α, so that (e^( -i k x))->(e^( -α x)); normalizability then requires B,C=0 but A,D≠0, which is only possible if S has a zero eigenvalue, rather than a pole, or equivalently  (S^(-1)) has a pole.  Does this give the same result?  Let’s check.  First, let’s check the condition that (S^(-1))=S†  has a pole:

Collect[Denominator[S\[ConjugateTranspose][[1,1]]//ComplexExpand//Simplify]/.{η->-I β}//TrigToExp,{(g^2)}]==0

This is precisely what we got above.  Now let’s check the condition that S has a zero eigenvalue.  A simple test for this is that the determinant must vanish.

Numerator[Det[S/.{η->-I β}//TrigToExp]//Simplify]==0

Which, again, is the same condition.  Good.  The lesson here is that we could use either sign in the analytic continuation η->± i β, but which sign we choose determines which of {A,D} or {B,C} must vanish for normalizability, and thus whether S must have a pole or a zero.