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

Question

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:

Answer 1

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

Answer 2

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'
        \addplot [red,samples=49] {V(x)};

        \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} {
            % add the "wave" lines
            \addplot [
                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'
            \addplot [
                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}