Intersection Between Stable and Saddle-Point Solutions -Stability Analysis

Question

I have a free energy function:

$$G(N_b, l_b)= -N_b E_b + \frac{1}{2} N_b \kappa_b (l_b - 1)^2 - F(l_b - 1) + \frac{A}{2} k_g (u_g - l_b)^2 + (N_t - N_b) \text{Log}\big(\frac{N_t - N_b}{A}\big) + N_b \text{Log}\big(\frac{N_b}{A}\big)$$

Physically, this represents the fee energy function of two substrates conjugated by $N_b$ reversible bonds that are being pulled away by a constant force $F$. Minimizing the free energy $G$ with respect to the state variables $N_b$ and $l_b$ ($\partial G/\partial N_b=0$ and $\partial G/\partial l_b=0$) yields the equilibrium equations for the system. I have solved these equations numerically. For a given force $F$, I obtain two roots, one is the stable solution (solid-line branch) and the other one is the saddle-point solution (dashed-line branch) as shown in the figure below. I used the second derivative test (checking the signs of Hessian eigenvalues) to make sure these roots are indeed the local minima and saddle-point in the free energy landscape.

You can see that at a critical force $F_{cr}$ these two branches intersect, this is the force at which the assembly snaps and the two substrates completely separate from each other. What I want to obtain is an (analytical) expression for $F_{cr}$ as a function of one of the input variables, such as $\kappa_b$. More specifically, my question is that mathematically what happens at the intersection. Of course, finding an analytical expression of $F_{cr}$ as a function of $\kappa_b$ would be ideal but using a numerical method to get $F_{cr}$ vs $\kappa_b$ curve would be great too. The curves shown above are constructed using $N_t=50$, $E_b=10$, $k_g=2.43902\times10^{-6}$, $u_g=3$, and $A=100$.

Answer 1

With

$$ G(N_b,l_b,F,\kappa_b) $$

calculating the partial derivatives

$$ \cases{ G_{N_b} = g_1(N_b,l_b,F,\kappa_b)\\ G_{l_b} = g_2(N_b,l_b,F,\kappa_b) } $$

and extracting $N_b$ from $g_2 = 0$ and substituting int $g_1=0$ we obtain

$$ \hat g_1(l_b,F,\kappa_b) = 0 $$

Plotting this function for $\kappa_b = 0.5,1.5,\cdots,9.5$ we obtain a set of blue curves as depicted in the following plot.

To obtain the red curve representing $F_{cr}(F,l_b)=0$ we proceed as follows:

$$ \frac{\partial \hat g_1}{\partial N_b}dN_b+\frac{\partial \hat g_1}{\partial l_b}dl_b=0 $$

now

$$ \frac{dN_b}{d l_b} = -\frac{\frac{\partial \hat g_1}{\partial l_b}}{\frac{\partial \hat g_1}{\partial N_b}} = g_3(l_b,F,\kappa_b)=0 $$

Extracting $\kappa_b$ from $g_3=0$ and substituting into $\hat g_1(l_b,F,\kappa_b)$ we obtain $\hat{\hat{g}}_1(l_b,F_{cr})=0$. This last function, generates the red dashed curve previously shown.

Attached a MATHEMATICA script to perform those transformations.

parms = {nt -> 50, eb -> 10, kg -> 2.43902 10^-6, ug -> 3, a -> 100};
g = -nb eb + 1/2 nb kb (lb - 1)^2 - f (lb - 1) + a/2 kg (ug - lb)^2 + (nt - nb) Log[(nt - nb)/a] + nb Log[nb/a];
equs = Grad[g, {nb, lb}];
solnb = Solve[equs[[2]] == 0, nb][[1]];
equs1 = equs[[1]] /. solnb;
grad = Grad[equs1, {lb, f}];
dlbdf = grad[[1]]/grad[[2]] /. parms // FullSimplify // Numerator // Chop;
solkb = Quiet@Solve[dlbdf == 0, kb];
equs11 = equs1 /. solkb[[2]] /. parms;
grf2 = Quiet@ContourPlot[equs11 == 0, {f, 10, 600}, {lb, 1, 6}, ContourStyle -> {Red, Dashed}]
equs10 = equs1 /. parms;
grf0 = Show[Table[ContourPlot[equs10 == 0, {f, 10, 600}, {lb, 1, 6}, PlotPoints -> 200], {kb, 0.5, 10}]];
Show[grf0, grf2]