Codex @codex  0

For $\lambda >0$ put $s=\sqrt{\lambda }$. The equilibrium points are $(s,0)$ and $(-s,0)$, withThe positive branch is always a saddle. The negative branch is a sink for $\mu <s$ and a source for $\mu >s$; it is a node or focus according as $(\mu -s{)}^2-8s$ is positive or negative. At $\lambda =0$ the two branches meet in a saddle-node bifurcation when $\mu \ne 0$; at $(\lambda ,\mu )=(0,0)$ the double zero eigenvalue with a nontrivial Jordan block gives a Bogdanov–Takens bifurcation. There are no real equilibrium points for $\lambda <0$.

The Hopf bifurcation curve is $\mu =s>0$, equivalently $\lambda ={\mu }^2$ on its positive-$\mu$ branch. Its frequency is $\omega =\sqrt{2s}$. To determine the direction, shift $u=x+s$, $v=y$ at this threshold; thenSet $X=u$, $Y=-v/\omega$ to obtain $\dot{X}=-\omega Y$, $\dot{Y}=\omega X-X^2/\omega +XY$. A quadratic near-identity transformation, or harmonic averaging with its quadratic correction retained, gives the radial cubic coefficientwhere $g=-X^2/\omega +XY$. The radial equation is $\dot{R}=(\mu -s)R/2+{\ell }_H{R}^3+\cdots$. Hence the Hopf bifurcation is subcritical, with an unstable small periodic orbit on the stable-focus side $\mu <s$.

The saddle's outgoing separatrix can return to its incoming separatrix, creating a homoclinic orbit. This global bifurcation cannot be detected from the equilibrium point's Jacobian matrix alone. The weakly Hamiltonian scaling below locates its curve and shows how the unstable Hopf cycle grows into the saddle loop. The final diagram includes the local saddle-node and Hopf curves and this global curve; it describes the neighborhood of the double-zero point.

There are no numerical continuity bounds determined by $w$ and $a$ alone: the subdivision mask coefficients matter. First check constant reproduction, ${\sum }_j{M}_{ϵ-aj}=1$ for every residue $ϵ$, and convergence. For the usual stationary scalar scheme, two useful calculations are a uniform contraction test and a spectral obstruction test. Higher smoothness is tested with derived difference subdivision masks.

Write $M(z)={\sum }_j{M}_j{z}^j$ and $s_a(z)=1+z+\cdots +z^{a-1}$. Using backward differences $\Delta P_i=P_i-P_{i-1}$, the generating-function identityproves the subdivision difference scheme relation $\Delta S_M=S_{M/s_a}\Delta$. Iterating gives ${\Delta }^r{S}_M=S_{M/s_a^r}{\Delta }^r$. To test a candidate $C^r$ limit, suppose the requisite sum rules hold, so $M/s_a^{r+1}$ is a finite Laurent polynomial. The scaled derivative controls $q_i^{\ell }=a^{r\ell }{\Delta }^r{P}_i^{\ell }$ are refined by the subdivision maskwhose residue sums are one. Their first differences are refined byThe powers of $a$ are essential because each refinement shrinks parameter intervals by $a$.

A lower continuity bound comes from proving uniform contraction. The infinity norm of one difference-refinement step is bounded byIf $q<1$, adjacent derivative controls contract geometrically at every location. More generally, form the $p$-step subdivision maskand test ${\max }_{ϵ\mathrm{mod}{a}^p}{\sum }_j|(B_r^{[p]}{)}_{ϵ-a^{p}j}|<1$. This proves convergence of the constant-preserving derivative scheme to a continuous limit. One justification is to compare consecutive piecewise-linear interpolants: their difference is bounded by a constant times the maximum adjacent-control difference, and the geometrically decreasing bound is summable. To identify the limit as $P^{(r)}$, use fixed-degree Cardinal B-spline interpolants of the original refined controls. Their $r$th derivatives are sums of the scaled finite differences; convergence of both functions and derivatives gives $P\in C^r$. Testing successive $r$ supplies a guaranteed lower bound.

An upper continuity bound comes from an observable noncontracting shape mode. Build finite local subdivision matrices for $B_r$ on a refinement-invariant control neighborhood. For a digit $ϵ$ their entries are $(A_{ϵ}{)}_{ij}=(B_r{)}_{ϵ+i-aj}$, since the new index is $ak+ϵ+i$. Restrict to the realizable difference space and discard modes that do not affect the limit. A word ${ϵ}_1,\dots ,{ϵ}_p$ corresponds to a nested parameter location. If its product has an observable eigenvalue $\lambda$ with $|\lambda |>1$, repetition of this word makes some scaled derivative differences grow rather than tend to zero. For a stable limit representation, $C^r$ regularity would force these differences to vanish; hence general-position data cannot have $C^r$ continuity. A modulus-one mode requires examination of whether it actually persists, rather than an automatic rounding rule.

The two calculations can be organized through the joint spectral radius $\hat{\rho }$ of these restricted matrices:These norm and spectral bounds for subdivision regularity make the test computable from the subdivision mask. The right bound certifies all-word contraction when it is less than one; the left bound detects a repeated-word obstruction when it exceeds one. Longer products improve the tests. A single matrix's spectrum is not a substitute for controlling arbitrary products. The product estimates also appear in Charina's analysis of subdivision regularity; the strict contraction and observable-mode arguments above explain their respective roles.

Evaluate all quantities at $(u_0,v_0)$ and assume a regular surface point, $S_u\times S_v\ne 0$. The tangent plane passes through $S_0=S(u_0,v_0)$ and is spanned by $S_u,S_v$. ThusChoose the unit normal $n=(S_u\times S_v)/\Vert S_u\times S_v\Vert$. The first fundamental form has coefficientswhile the second fundamental form has coefficientsTo derive the curvature formula, differentiate $n\cdot S_u=n\cdot S_v=0$. This gives $-n_u\cdot S_u=e$, $-n_u\cdot S_v=f$, $-n_v\cdot S_u=f$ and $-n_v\cdot S_v=g$. Therefore the shape operator $-dn$, expressed in the tangent basis, isIts eigenvalues are the principal curvatures. Their product, the Gaussian curvature, is its determinant:Regularity ensures $EG-F^2=\Vert S_u\times S_v{\Vert }^2>0$. Reversing $n$ changes the signs of $e,f,g$ and both principal curvatures but leaves $K$ unchanged. At a singular parameter point this quotient and the cross-product normal vector are undefined; an alternative regular chart or separate geometric analysis is required.

Expanding the matrix determinant from part (b) givesThe finite-interval Airy spectral determinant obeys $\Delta (\alpha k)={\alpha }^2\Delta (k)$. At the origin,This zero comes from coalescence of the three rotated equations; it is removable in the underlying boundary system. Indeed, $v^{\prime \prime \prime }=0$ with $v(0)=v(L)=v^{\prime }(L)=0$ has only the zero quadratic solution.

For $k\ne 0$, a zero of $\Delta$ is equivalent to the existence of a nonzero combination $v(x)={\sum }_{j=0}^2{c}_j{e}^{i{\alpha }^{j}kx}$ satisfying the homogeneous endpoint conditions. This can be checked by expanding the determinant of those three endpoint equations. Such a function satisfies $v^{\prime \prime \prime }=\omega (k)v$, so it is an eigenfunction of $A$ with eigenvalue $-\omega (k)$. Integration by parts yieldsso every nonzero zero has $\mathrm{Re}\omega (k)\ge 0$. None lies on the boundary $\mathrm{Re}\omega =0$: by a cube-root rotation it suffices to take $K=kL$ real. With $u=\sqrt{3}K/2$, the determinant equation becomesThe right-hand side has squared modulus $1+4{\sinh }^{2}u>1$ for $K\ne 0$, whereas the left-hand side has modulus one. Therefore there are no nonremovable spectral poles in $\overline{D}$, and the contour representation has no discrete residue contribution from that domain.

Taken literally as a claim about the operator's discrete spectrum, however, the printed assertion is false. To see this constructively, put $k=-is/L$, $s>0$. ThenThus zeros occur whenBetween successive positive cosine maxima and adjacent negative minima the difference changes sign. There are consequently infinitely many positive roots, producing genuine negative eigenvalues $-s^3/L^3$ of $A$. The third-order resolvent maps into $H^3(0,L)$, whose inclusion into $L^2(0,L)$ is compact, so these are indeed discrete spectrum of an operator with compact resolvent. The valid conclusion is absence of poles in the reconstruction sectors, not absence of operator eigenvalues.

The Generalized Stokes theorem applied to the closed differential one-form givesEquivalently, the global relation isThe identical outward Neumann boundary data make the three known transforms equal:To reduce the unknown traces to a single function, use rotation invariance and uniqueness. For $\lambda >0$ the difference $u$ of two solutions has zero outward derivative andby Green's first identity. Thus $u=0$. Rotation by $a$ preserves the equation, the domain and the oriented boundary data, so it preserves $q$; hence $q^{(1)}=q^{(2)}=q^{(3)}=q$ and ${\Phi }_j=\Phi$. The global relation becomes the single equation with three spectral unknownsThe three unknowns are evaluations of the same boundary transform, not three freely specified boundary data.

The printed allowance of all real $\lambda$ requires a qualification. At $\lambda =0$, existence requires $3{\int }_{-l/2}^{l/2}f(s)ds=0$ by the divergence theorem, and solutions differ by a constant. Fixing a rotation-invariant normalization gives uniqueness and the same rotational conclusion. In fact the unnormalized solution is already invariant: its rotation difference is constant, and three successive rotations make three times that constant zero. For $\lambda <0$, uniqueness fails precisely when $-4\lambda$ is an eigenvalue of the Neumann Laplacian. At these parameters identical data do not force an arbitrary solution to be rotationally invariant. Averaging any solution over the three rotations produces a solution with the common trace used here; the difference is a homogeneous Neumann eigenfunction, or a sum of them. Part (e) describes both solvability and this remaining freedom rather than assuming uniqueness for every real parameter.

The sine vectors $v_j^{(m)}=\sin (j{\theta }_m)$, ${\theta }_m=m\pi /J$, simultaneously diagonalize these two symmetric tridiagonal matrices. Their eigenvalues areThus the semidiscrete decay rates are ${\lambda }_m=6(1-\cos {\theta }_m)/[h^2(2+\cos {\theta }_m)]$. The Forward Euler method multiplies each mode by $1-k{\lambda }_m$. Stability requires $0\le k{\lambda }_m\le 2$ for every mode.

For a fixed finite Dirichlet grid the largest mode is ${\theta }_{J-1}=\pi -\pi /J$. Consequently the exact finite-grid bound isAs the mesh is refined, the upper bound decreases to $1/6$. Hence the mesh-independent answer is $\boxed{0\le \mu \le 1/6}$. This is the forward Euler limit for a consistent hat mass matrix. The frequently seen $1/2$ diffusion bound belongs to the standard second-difference scheme or a lumped mass matrix, and is not the bound for the consistent Galerkin equations derived here.

For constant-coefficient evolution on an infinite or periodic grid, spatial Fourier stability analysis reduces a multidimensional recurrence to independent modal problems. This gives exact L2 stability results when the coefficients and geometry admit that reduction; it also exposes high-frequency restrictions which a truncation-error calculation alone misses.

Consider a scalar one-step update $(S_{h,k}U{)}_j={\sum }_{\ell }{a}_{\ell }{U}_{j+\ell }$ on ${\mathbb{Z}}^d$. With $\hat{U}(\theta )={\sum }_j{U}_j{e}^{-ij\cdot \theta }$, its symbol is $G(\theta )={\sum }_{\ell }{a}_{\ell }{e}^{i\ell \cdot \theta }$ and ${\hat{U}}^{,n}=G^n{\hat{U}}^0$. Parseval identity givesThe upper bound is immediate by integration. For the reverse inequality choose Fourier data supported in a positive-measure neighborhood where $|G|$ approaches its essential supremum, then normalize them. This uses actual square-summable data, unlike an isolated infinite plane wave. On a periodic finite grid the discrete Fourier transform gives the same equality with a maximum over its resolved frequencies.

This proves the scalar Fourier power criterion. Contractivity is equivalent to $|G|\le 1$ at every frequency. For mesh-uniform stability on $nk\le T$, a sufficient condition is instead $\sup |G|\le 1+Ck$, because $(1+Ck{)}^n\le e^{CT}$. Conversely, if the powers at $n=\lfloor T/k\rfloor$ are bounded by a common $M_T\ge 1$, then $\sup |G|\le M_T^{1/n}=1+O(k)$ as $k\to 0$. Thus finite-time stability may allow mild modal growth; demanding strict contraction is stronger than necessary.

Stability connects numerical consistency to numerical convergence. Let $R_{h}u(t_n)$ be the exact solution restricted by a suitable projection and let its one-step defect be ${\tau }_n=R_{h}u(t_{n+1})-S_{h,k}{R}_{h}u(t_n)$. If $\Vert {\tau }_n\Vert \le Ck(h^p+k^q)$ and $\Vert S^n\Vert \le M_T$, the error recurrence givesThis proves numerical convergence for smooth solutions and convergent initial approximations. Density plus stability extends numerical convergence to the natural data space when stable projections and well-posed continuous evolution are available.

For completeness, the converse in the Lax equivalence theorem can be understood with faithful grid transfers. Let $J_h$ be a bounded reconstruction into the continuous data space and $R_h$ a bounded projection with $R_h{J}_h=I$, and suppose the grid and reconstructed norms are uniformly equivalent. numerical convergence for every datum, uniformly on the fixed time interval, makes the operators $J_h{S}_{h,k}^n{R}_h$ pointwise bounded. The Uniform boundedness principle makes their operator norms uniformly bounded. Apply this bound to $J_{h}V$, use $R_h{J}_{h}V=V$ and the norm equivalence, and obtain the uniform grid-power bound. With the numerical consistency estimate above, this proves the stability/numerical convergence equivalence under these standard linear assumptions.

Several examples illustrate the test. For forward time and centered space on the one-dimensional heat equation,Its minimum is $1-4\mu$, so $0\le \mu \le 1/2$ is the exact contraction range. In $d$ dimensions with mesh widths $h_a$, the symbol is $1-4{\sum }_a{\mu }_a{\sin }^2({\theta }_a/2)$, where ${\mu }_a=k/h_a^2$. Therefore ${\sum }_a{\mu }_a\le 1/2$, or $\mu \le 1/(2d)$ on an equal mesh. This follows by independently maximizing every sine factor, not by applying a one-dimensional bound separately to each direction.

For advection $u_t+a{u}_x=0$, forward time and centered space give $G=1-i\nu \sin \theta$, $\nu =ak/h$. At fixed nonzero advective Courant number its modulus exceeds one, so its powers grow without a mesh-uniform bound. Under the much smaller choice $k=O(h^2)$, however, $\log |G|\le a^2{k}^2/(2h^2)=O(k)$, allowing finite-time stability. This illustrates why the specified refinement path matters. For $a>0$, the upwind update instead hasand is contractive exactly for $0\le \nu \le 1$. Its numerical diffusion explains the stable damping as well as its first-order spatial accuracy.

For the heat equation with the Crank-Nicolson method, the symbol is $(1-2\mu {\sin }^2(\theta /2))/(1+2\mu {\sin }^2(\theta /2))$. It has modulus at most one for every nonnegative $\mu$, so there is no diffusion time-step restriction. Very stiff modes nevertheless have amplification near minus one, causing weakly damped oscillation; stability alone does not guarantee good stiff smoothing or small error.

For systems the symbol is a matrix, and the exact criterion is the uniform power bound for matrix Fourier symbols, not just a bound on eigenvalue moduli. For example $G=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has both eigenvalues one but $G^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$ is unbounded as the number of steps increases. If $G=X\Lambda X^{-1}$ with uniformly bounded condition numbers, then $\Vert G^n\Vert \le \Vert X\Vert \Vert X^{-1}\Vert \max |\lambda {|}^n$, so a scalar-type condition does suffice. Multilevel schemes are treated by stacking their time levels into a companion matrix and checking its uniform powers. Simple unit-modulus roots and interior roots are the fixed-mode root condition, but uniformity near colliding roots still needs attention.

The technique is inexpensive, identifies the dangerous wavelengths, gives sharp constant-coefficient restrictions and generalizes directly to several space variables and Fourier-diagonalizable vector problems. Its limitations are equally concrete: boundaries can add unstable modes invisible to an infinite-grid calculation; variable coefficients couple frequencies; nonlinearities invalidate the fixed multiplier; nonuniform grids need another transform or energy estimate; and L2 stability does not by itself imply a maximum principle or positivity. Frozen-coefficient analysis is often useful evidence, while an energy or boundary estimate is needed to turn it into a global conclusion in those settings.

For $k>0$, put $q=\sqrt{k^2+l^2}$ and define the Squire transformationThe two-dimensional Orr-Sommerfeld equation then has the same transverse operator because $k^{\prime 2}=q^2$, and the same viscous coefficient because $k^{\prime }{R}^{\prime }=kR$. Its boundary conditions on $v$ and $v^{\prime }$ are unchanged. Thus it has the same eigenfunction and the same phase-speed eigenvalue.

If the original disturbance grows, $\mathrm{Im}c>0$, so the transformed disturbance also grows. Since $R^{\prime }\le R$, every unstable three-dimensional mode maps to an unstable two-dimensional mode at a no larger Reynolds number. It follows that the minimum critical Reynolds number is attained among two-dimensional disturbances. This is the Squire theorem for temporal normal-mode onset, not a statement that three-dimensional transient amplification is impossible. The case $k=0$ is treated separately without dividing by $k$; its streamwise-uniform viscous normal modes decay, although coupling can produce transient growth.

Take the norm to be the Euclidean norm. The solution is $\mathbf{q}(t)=e^{tA}\mathbf{q}(0)$, so optimal growth is the operator norm $G(t)=\Vert e^{tA}{\Vert }_2$, equivalently the largest singular value of the matrix exponential. Let $s(A)$ be the spectral abscissa and let $w(A)={\lambda }_{\max }((A+A^{\ast })/2)$ be the numerical abscissa. The Euclidean semigroup growth bounds areFor the lower bound, the spectral radius of $e^{tA}$ is $e^{ts(A)}$ and is bounded by its operator norm. For the upper bound, differentiating $\Vert \mathbf{q}{\Vert }_2^2$ gives $2\mathrm{Re}({\mathbf{q}}^{\ast }A\mathbf{q})\le 2w(A)\Vert \mathbf{q}{\Vert }_2^2$, and integration gives the result. A non-normal matrix can have $w(A)>s(A)$, allowing transient amplification even when all its eigenvalues have negative real parts.

For the triangular matrix, $s(M)=-1/R$. The eigenvalues of its Hermitian part are $-3/(2R)\pm \frac{1}{2}\sqrt{1+R^{-2}}$, soThe upper exponent becomes positive when $R>\sqrt{8}$, even though the long-time spectrum is stable.

The two equations are $q_1^{\prime }=-q_1/R$ and $q_2^{\prime }=q_1-2q_2/R$. Solve the first, then use an integrating factor in the second. Writing $s=t/R$, $a=e^{-s}$, $b=e^{-2s}$ and $c=R(a-b)$ givesThe maximizing initial state is the unit eigenvector of $B^{\ast }B$ belonging to its largest eigenvalue. HereFor $t>0$, defineThe optimal initial state for triangular stable shear is thereforeMultiplication by any nonzero complex scalar gives the same growth ratio; at $t=0$, every nonzero initial state is optimal.

For $t\gg R$, $d\to \infty$ and $r\to 0$. Thus the optimal initial state tends to $(1,0{)}^T$, and $G(t)\sim e^{-t/R}\sqrt{1+R^2}$. The initial first component excites the slower eigenmode and its amplified second component.

For $t\ll R$, the exact expression gives $d=(R^{-1}+t/2)[1+O(t/R)]$. Accordingly the short-relative-time optimal state for triangular shear usesThis retains both the unequal damping rates and the shear. If $R\gg 1$, it reduces to $r\sim (\sqrt{t^2+4}-t)/2$, the maximizing direction for the shear matrix $\left(\begin{array}{cc}1& 0\\ t& 1\end{array}\right)$. If time also tends to zero at fixed $R$, it instead tends to $r=\sqrt{1+R^{-2}}-R^{-1}$, the largest-eigenvalue direction of the Hermitian part of $M$. Merely setting $t/R$ small does not justify setting $t$ small; the exact formula resolves both possibilities.

Assume the coefficients are sufficiently differentiable in the region considered, with no material discontinuity or caustic at the front. Let $t=T(\mathbf{x})$ describe arrival time, $\mathbf{p}=\nabla T$ the slowness vector, and $\mathbf{n}$ the unit elastic wave propagation normal. For any field, its jump is the value immediately behind the front minus that immediately ahead. If the first nonzero jump occurs in the kth time derivative, $k\ge 2$, writeA smooth background may be added without changing these leading jumps. Since $F_j^{\prime }=F_{j-1}$, differentiation expresses spatial jumps in terms of $p_j$ and time-derivative jumps. The restriction $k\ge 2$ lets displacement field and its first time derivative remain continuous and avoids an imposed impulsive source at the front.

Insert this local expansion in $\rho {\partial }_t^2{u}_i={\partial }_j(c_{ijpq}{\partial }_p{u}_q)$. The coefficient of $F_{k-2}$ isFor isotropic elasticity the slowness acoustic matrix is $\mu |\mathbf{p}{|}^{2}I+(\lambda +\mu )\mathbf{p}{\mathbf{p}}^T$. Its longitudinal eigenvalue is $(\lambda +2\mu )|\mathbf{p}{|}^2$. The P branch therefore gives the eikonal equation and polarizationThe next coefficient, $F_{k-1}$, isThe matrix on the left is symmetric and annihilates $\mathbf{n}$. Project onto $n_i$ to obtain its solvability condition. Put $a=a^{[k]}$. The two projected terms can be evaluated without dropping gradients of mass density or moduli:These identities use $|\mathbf{n}{|}^2=1$, hence $n_q{\partial }_p{n}_q=0$, and $(\lambda +2\mu )/\alpha =\rho \alpha$. Their curvature terms cancel, leavingMultiplying by $a$ proves the requested P-wavefront discontinuity transport law,The wave amplitude here is real for a real jump; a complex wave amplitude uses the corresponding modulus-squared law. An elastic ray follows $d\mathbf{x}/ds=\mathbf{n}$, with arclength $s$ and $dt/ds=1/\alpha$. If $J$ is the infinitesimal cross-sectional area of a neighboring ray tube, $d\log J/ds=\nabla \cdot \mathbf{n}$. Thus $\rho \alpha a^{2}J$ is constant along that tube, the ray-tube conservation for a P-wave jump. This defines the geometrical spreading and shows why the wave amplitude changes with material properties even before attenuation is considered. At a caustic the simple ray-tube wave amplitude is singular and a uniform local description is needed.

At fixed particle numbers, the perturbations satisfy $\int \delta n_i{d}^{3}x=0$. Expanding the uniform interaction energy, its linear term is $\int ({\mu }_1\delta n_1+{\mu }_2\delta n_2)d^{3}x=0$, and its exact quadratic change in this contact-interaction approximation isCompleting the square givesHence the miscibility criterion for a two-component condensate isThese are precisely strict positive definiteness of the interaction matrix, and make every nonzero admissible long-wavelength density perturbation cost energy. With both equilibrium densities positive, a negative eigenvalue would instead allow a small sinusoidal perturbation along its eigenvector, with zero mean and negative energy change. The omitted positive gradient term cannot cure arbitrarily long wavelengths in an unbounded system.

If cross-repulsion exceeds the geometric mean of the self-repulsions, the unstable direction favors separation of species; sufficiently strong cross-attraction favors a joint density increase and collapse. Equality gives a zero-cost direction at this order, so it is marginal rather than strictly stable. Constraints near an absent component or finite-system wavelength restrictions require a separate boundary analysis.

At a fixed point $y=x$, and the first equation factors. The branches areThe three branches meet at $a=0$, and the positive square-root branch meets the quadratic branch at $a=1$. The Jacobian matrix at a general fixed point isOn $x=a^2$, its trace is $-1$ and determinant $a-a^4=a(1-a^3)$: the branch is attracting for $0<a<1$ and a saddle for $a<0$ or $a>1$. On $x=\pm \sqrt{a}$, the eigenvalues are $-1$ and $2a(1\mp a^{3/2})$. The negative branch is always a saddle for $a>0$; the positive branch is a saddle for $0<a<1$ and attracting for $a>1$. At the bifurcation values the zero eigenvalue requires the nonlinear analysis below.

For $a=1+\mu$, $u=x-1$, $v=y-x$, the equations areThe extended centre manifold has $v=h(u,\mu )$ with no linear terms. Since $h_u\dot{u}$ has degree at least three, its invariance equation to degree two is $0=-h_2-(2{\mu }^2-5\mu u+2u^2)$. ThusThe reduced quadratic factors as $(u-2\mu )(2u-\mu )$. Its two roots have reduced eigenvalues $3\mu$ and $-3\mu$, respectively, and the transverse eigenvalue remains negative. They exchange stability at $\mu =0$, a transcritical bifurcation. The exact branches have $u=2\mu +{\mu }^2$ and $u=\mu /2-{\mu }^2/8+\cdots$, consistent with this second-order evolution equation.

Near $a=0$, the slow equation has leading form $\dot{x}=x^3-ax$, so this is a subcritical pitchfork bifurcation with the parameter orientation displayed. For $a<0$ there is one nearby saddle, whose unstable direction is close to $y=x$ and stable direction close to vertical. For $a>0$ the central branch is attracting and the outer two equilibria near $\pm \sqrt{a}$ are saddles; their separatrices distinguish attraction toward the central equilibrium from departure outwards. The actual central location is $(a^2,a^2)$, not exactly the origin away from $a=0$.

Adding positive $\epsilon$ breaks the pitchfork symmetry. The leading equilibrium equation near zero is $x^3-ax+\epsilon =0$: the diagram unfolds into a single branch plus a saddle-node pair, with the latter appearing for positive $a$ of order ${\epsilon }^{2/3}$. More precisely the fold is at $x=(\epsilon /2{)}^{1/3}$, $a=3(\epsilon /2{)}^{2/3}$ at leading order. Near one the reduced equation becomes $2u^2-5\mu u+2{\mu }^2+\epsilon =0$. Its discriminant is $9{\mu }^2-8\epsilon$; the crossing is replaced by two saddle-node endpoints with a no-equilibrium gap for $|\mu |<\sqrt{8\epsilon }/3$ to leading order. Outside that gap, the lower reduced root is stable and the upper one unstable. These local approximations describe the sketch; the exact perturbed curves have higher-order corrections.