# How to plot transition due to Franck-Condon principle using pgfplots?

There are several ways to plot Morse potential (see this question) and other smart solutions to plot the quantum harmonic oscillator wave functions on a parabola (see this other question).

So far I copied the solution for the Morse potential as it looked very nice (I prefer when the machine execute loop instead of me coping 100 wave functions by hand):

% ... %
\usepackage{tikz}
\usetikzlibrary{intersections}
\usepackage{pgfplots}
\usepgfplotslibrary{fillbetween}

% ... %

\begin{tikzpicture}
%%%%%%% Define Potential Function %%%%%%%
\pgfmathsetmacro{\De}{6}
\pgfmathsetmacro{\Ro}{1}
\pgfmathsetmacro{\alpha}{1}
\pgfmathdeclarefunction{V}{1}{%
\pgfmathparse{%
\De*((1-exp(-\alpha*(#1-\Ro)))^2-1)%
}%
}%
%%%%%%% Energy Levels %%%%%%%
\pgfmathdeclarefunction{energy}{1}{%
\pgfmathparse{%
-\De+(#1+.5) - (#1+.5)^2/(4*\De)
}%
}%

\begin{axis}[
axis lines=none,
smooth,
no markers,
domain=0:8,
ymax=6,
scale=1.5
]
\addplot [red, samples=50, name path global=MorseCurve] {V(x)};
\pgfplotsinvokeforeach{0,1,2,3,4,5,6, 7, 8, 9, 10, 11, 12}{
\path [name path global=HelperLine-#1] (axis cs: 0,{energy(#1)}) -- (axis cs: 10, {energy(#1)});
\draw[name intersections={of=MorseCurve and HelperLine-#1}] (intersection-1) -- (intersection-2);
}
\end{axis}
\end{tikzpicture}
% ... %

Since I would like to represent the Franck-Condon principle, what I'm looking for is a way to merge the two things, the quantum harmonic oscillator anche the morse potential, so that I could achieve something like

I really would like to keep everything in a for loop-like structure, but I understand that it's impossible to compute Hermite polynomials with LaTeX (isn't it?) so maybe a labelling with the numbers of \pgfplotsinvokeforeach would even be a great step forward. Any idea?

# EDIT 1

The first 4 functions are (as requested)

H_1(x) = exp(-x^2/2)
H_2(x) = 2*x * exp(-x^2/2)
H_3(x) = (4*x^2-2) * exp(-x^2/2)
H_4(x) = (-12*x+8*x^3) * exp(-x^2/2)
...

The polynomial part are the Hermite polynomials, while the exponential part it just remains the same. Here I computed them with Mathematica up to 10:

• I think there maybe a solution to your problem, but for that of course need either the function to calculate the waves or (even better) a data file or data files of them. Could you provide these? Commented Jun 3, 2018 at 16:30
• The function I guess you'll find impossible to reproduce on pgf because it's defined as the product of the hermite polynomials multiplied by a gaussian. So I'll provide the analytical form, and then I'll try produce a datafile. Commented Jun 3, 2018 at 17:43
• I think producing a data table/file when you have the function should be a simple task. I plotted H1 and H10 using gnuplot and noticed that the the max value of H1 is 1 and the max value of H10 is around 4*10^4. So when you produce the file you should directly somehow "normalize" the height of the waves to be drawn in a proper way. I would suggest a height of 1. If there need to be some additional scaling this can be done by PGFPlots. Commented Jun 3, 2018 at 18:13
• Then I think the waves are also shifted in x direction to the minimum value of the V function, right? Please also provide this value or directly shift the x values in the data table. Commented Jun 3, 2018 at 18:14
• Did one of the answers help you to solve the question? If yes, please consider upvoting (by clicking on the arrows next to the score) and/or marking it as the accepted answer (by clicking on the checkmark ✓). Commented Jun 10, 2018 at 19:21

It is quite hard, at least for me, to do that with pgfplots because of the various expansion issues. You only provide the first four of the Hermite polynomials, and do not normalize them. Therefore, I present a TikZ "only" solution in which I make the normalizations up. But I will be happy to improve this solution once you provide all relevant polynomials including normalization.

\documentclass[tikz,border=3.14mm]{standalone}
\usetikzlibrary{intersections,calc}
\begin{document}
\begin{tikzpicture}[yscale=0.5]
\tikzset{declare function={H0(\x)=pi^(-1/4)*exp(-\x*\x);
H1(\x)=(sqrt(2)*\x)/pi^(1/4)*exp(-\x*\x);
H2(\x)=(-(1/(sqrt(2)*pi^(1/4))) +  (sqrt(2)*abs(\x^2))/pi^(1/4))*exp(-\x*\x);
H3(\x)=(-((sqrt(3)*\x)/pi^(1/4)) +
(2*\x^3)/(sqrt(3)*pi^(1/4)))*exp(-\x*\x);
H4(\x)=(sqrt(3)/(2*sqrt(2)*pi^(1/4)) -
(sqrt(6)*abs(\x^2))/pi^(1/4) + (sqrt(2/3)*abs(\x^4))/pi^(1/4))*exp(-\x*\x);
H5(\x)=((sqrt(15)*\x)/(2*pi^(1/4)) -
(2*sqrt(5)*\x^3)/(sqrt(3)*pi^(1/4)) +
(2*\x^5)/(sqrt(15)*pi^(1/4)))*exp(-\x*\x);
H6(\x)=(-sqrt(5)/(4*pi^(1/4)) + (3*sqrt(5)*abs(\x^2))/
(2*pi^(1/4)) - (sqrt(5)*abs(\x^4))/pi^(1/4) +
(2*abs(\x^6))/(3*sqrt(5)*pi^(1/4)))*exp(-\x*\x);
H7(\x)=(-(sqrt(35)*\x)/(2*sqrt(2)*pi^(1/4)) +
(sqrt(35/2)*\x^3)/pi^(1/4) - (sqrt(14)*\x^5)/
(sqrt(5)*pi^(1/4)) + (2*sqrt(2/35)*\x^7)/(3*pi^(1/4)))*exp(-\x*\x);
H8(\x)=(sqrt(35)/(8*sqrt(2)*pi^(1/4)) -
(sqrt(35)*abs(\x^2))/(sqrt(2)*pi^(1/4)) +
(sqrt(35/2)*abs(\x^4))/pi^(1/4) - (2*sqrt(14)*abs(\x^6))/
(3*sqrt(5)*pi^(1/4)) + (sqrt(2/35)*abs(\x^8))/(3*pi^(1/4)))*exp(-\x*\x);
H9(\x)=((3*sqrt(35)*\x)/(8*pi^(1/4)) -
(sqrt(35)*\x^3)/pi^(1/4) + (3*sqrt(7)*\x^5)/
(sqrt(5)*pi^(1/4)) - (4*\x^7)/(sqrt(35)*pi^(1/4)) +
(2*\x^9)/(9*sqrt(35)*pi^(1/4)))*exp(-\x*\x);
H10(\x)=((-3*sqrt(7))/(16*pi^(1/4)) +
(15*sqrt(7)*abs(\x^2))/(8*pi^(1/4)) - (5*sqrt(7)*abs(\x^4))/
(2*pi^(1/4)) + (sqrt(7)*abs(\x^6))/pi^(1/4) -
abs(\x^8)/(sqrt(7)*pi^(1/4)) + (2*abs(\x^10))/(45*sqrt(7)*pi^(1/4)))*exp(-\x*\x);}
}

%%%%%%% Define Potential Function %%%%%%%
\pgfmathsetmacro{\De}{6}
\pgfmathsetmacro{\Ro}{1}
\pgfmathsetmacro{\alpha}{1}
\pgfmathdeclarefunction{V}{1}{%
\pgfmathparse{%
\De*((1-exp(-\alpha*(#1-\Ro)))^2-1)%
}%
}%
%%%%%%% Energy Levels %%%%%%%
\pgfmathdeclarefunction{energy}{1}{%
\pgfmathparse{%
-\De+(#1+.5) - (#1+.5)^2/(4*\De)
}%
}%

\draw[red, name path global=MorseCurve]
plot[samples=100,variable=\x,domain=0.25:8]
({\x},{V(\x)});
\foreach \X in {0,...,5}{
\path [name path global=HelperLine-\X] (0,{energy(\X)}) --
(10, {energy(\X)});
\path[name intersections={of=MorseCurve and HelperLine-\X}] (intersection-1) -- (intersection-2);
\draw let \p1=(intersection-1),\p2=(intersection-2) in
\pgfextra{\pgfmathsetmacro{\h}{\y1/1cm}\xdef\h{\h}
\pgfmathsetmacro{\xleft}{(\x1-1cm)/1cm}\xdef\xleft{\xleft}
\pgfmathsetmacro{\xright}{(\x2+1cm)/1cm}\xdef\xright{\xright}}
(\xleft,\h) -- (\xright,\h);
\pgfmathsetmacro{\xmid}{(\xleft+\xright)/2}
\filldraw[fill=orange] plot[domain=\xleft:\xright,variable=\x,samples=100]
({\x},{H\X(2.5*(\x-\xmid))/2+\h}) -- (\xright,\h) -- (\xleft,\h);
}
\end{tikzpicture}
\end{document}

UPDATE: Fixed an issue \x^2 vs. \x*\x, made the plots centered and added all normalization factors (using Mathematica). Unfortunately, I cannot plot the higher functions because of dimension too large errors'....

• Ok, I'll provide more function and normalize them, meanwhile thanks a lot! Commented Jun 4, 2018 at 11:38

I guess this is quite near what you want to achieve. Because a lot of stuff is done here, please have a look at the comments in the code to find out what is good for what.

(For fast compilation I recommend using LuaLaTeX.)

% used PGFPlots v1.16
\documentclass[border=5pt]{standalone}
\usepackage{pgfplots}
\usepgfplotslibrary{fillbetween}
\pgfplotsset{
% use this compat' level or higher to use Lua as calculation engine
% when compiling/TeXing with LuaLaTeX
compat=1.12,
}
\begin{document}
\begin{tikzpicture}[
% redeclared the variables/functions to a bit simpler form (as I think)
/pgf/declare function={
% Define Potential Function
De = 6;
Ro = 1;
alpha = 1;
V(\x) = De * ( (1 - exp(-alpha*(\x-Ro)))^2 - 1 );
%
% inverse function(s) of V(x)'
M(\V) = Ro - 1/alpha * ln( 1 + sqrt(\V/De +1) );
N(\V) = Ro - 1/alpha * ln( 1 - sqrt(\V/De +1) );
%        xVmin = Ro;         % which is the same as: M(-De); or N(-De);
%
% Energy Levels
energy(\z) = -De + (\z+0.5) - (\z+0.5)^2 / (4*De);
%
% how much should the horizontal lines be enlarged (in axis units)?
xshift = 1;
% ---------------------------------------------------------------------
% number of samples that should be used to plot the function
Samples(\n) = (1+\n)*5;
% common factor for Hn'
EE(\x)  = exp(-\x^2/2);
% (just "copied" from the image in the question)
H0(\x)  = EE(\x);
H1(\x)  = EE(\x) *   2*\x;
H2(\x)  = EE(\x) * (-2 + 4*\x^2);
H3(\x)  = EE(\x) * (-12*\x + 8*\x^3);
H4(\x)  = EE(\x) * (+12 - 48*\x^2 + 16*\x^4);
H5(\x)  = EE(\x) * (+120*\x - 160*\x^3 + 32*\x^5);
H6(\x)  = EE(\x) * (-120 + 720*\x^2 - 480*\x^4 + 64*\x^6);
H7(\x)  = EE(\x) * (-1680*\x + 3360*\x^3 - 1344*\x^5 + 128*\x^7);
H8(\x)  = EE(\x) * (+1680 - 13440*\x^2 + 13440*\x^4 - 3584*\x^6 + 256*\x^8);
H9(\x)  = EE(\x) * (+30240*\x - 80640*\x^3 + 48384*\x^5 - 9216*\x^7 + 512*\x^9);
H10(\x) = EE(\x) * (-30240 + 302400*\x^2 - 403200*\x^4 + 161280*\x^6 - 23040*\x^8 + 1024*\x^10);
%
% scaling factors to "normalize" the Hn' functions to around 1 (= max(Hn) - min(Hn))
% ("guessed" by manual inspection of the resulting graph/plot)
G0  = 1;
G1  = 2.4;
G2  = 4.3;
G3  = 4.6+5.4;
G4  = 12.5+15;
G5  = 38+46;
G6  = 130+160;
G7  = 480+600;
G8  = 1900+2300;
G9  = (0.8+0.95)*1e4;
G10 = (3.4+4.4)*1e4;
% additional scaling factor to further "shrink" the height
F  = 1.5;
},
]
\begin{axis}[
smooth,
domain=0:8,
axis lines=none,
ymax=6,
%        % ---------------------------------------------------------------------
%        % for debugging purposes only
%        % (comment the previous two lines and uncomment the following lines to
%        %  determine Gn')
%        minor y tick num=4,
%        xminorgrids,
%        % ---------------------------------------------------------------------
]
% adapted samples' so one of the sample points is near xVmin'

\pgfplotsinvokeforeach {0,...,10} {
% draw the horizontal lines at the different energy levels
% enlarge the lines by xshift'
\addplot [very thin,name path=h#1] coordinates {
({M(energy(#1)) - xshift},{energy(#1)})
({N(energy(#1)) + xshift},{energy(#1)})
}
% add the nodes to the right of the lines
node (m#1) [
node font=\footnotesize,
at end,
right,
xshift=6ex,
] {$#1$}
;
}

\pgfplotsinvokeforeach {0,...,6} {
name path=H#1,
very thin,
% adapt the domain' to consider the shift in Ro' and xshift'
% for the x values
domain={-Ro + M(energy(#1)) - xshift}:{-Ro + N(energy(#1)) + xshift},
% use a dynamic sample rate depending on the Hn'
samples={Samples(#1)},
] (
% shift x' by Ro'
{x + Ro},
% scale Hn' by Gn' and F', and shift the result by energy(n)'
{H#1(x)/G#1/F + energy(#1)}
);

% color the area between hn' and Hn'
fill=orange,
] fill between [
of=h#1 and H#1,
];
}

% add the missing label part of the m0'
\node [
node font=\footnotesize,
anchor=base east,
] at (m0.base east) {$v' = \phantom{0}$};

\end{axis}
\end{tikzpicture}
\end{document}

• I guess that in H9 an \x is missing. The Hn with odd n are odd functions....
– user121799
Commented Jun 4, 2018 at 20:42
• @marmot, of course, that makes perfect sense. Thank you. Added the missing \x to the first summand in the function definition. But I decided to only draw the "waves" up to H6 because then I can "enlarge" them a bit more. But users are free to choose what they like and only have to adjust F so the "waves" don't overlap any more. Commented Jun 4, 2018 at 21:14
• Your approach is better since I run into dimensions too large errors with mine. But I included the normalization factors of all functions in my post if you are interested....
– user121799
Commented Jun 4, 2018 at 21:23
• @marmot, thank you for the credits. I tested your "normalized" functions, but they seem to get flatter and flatter for greater |x| values, while OPs picture suggests that the outermost (shown) waves should be the highest. Therefore I think something is wrong with that normalization. Couldn't you just find the min and max values of the plotted region in Mathematica and scale by the difference of these two values? Commented Jun 5, 2018 at 5:51
• I used the quantum mechanics normalization since after all this is tagged physics. The quantum mechanics normalization says that \int dx |H_n(x)|^2=1. (The OP's H_n are not the Hermite polynomials, but the OP multiplies them by exp(-x^2/2). So I guess you may want to check with the OP what her or his preferences are. If (s)he wants a different normalization, I'll be happy to run Mathematica again, but I have a feeling (s)he won't want a different one.
– user121799
Commented Jun 5, 2018 at 16:22