The polar molecular field isFor a spatially uniform polar order parameter, $\mathbf{h}=(a+b{p}^2)\mathbf{p}$, so it is colinear with $\mathbf{p}$ and its magnitude is $|\mathbf{h}|=|a+b{p}^2|p$. The scalar factor can be negative; colinearity does not necessarily mean parallel orientation.

For the simple shear flow, the paper's velocity gradient convention givesSubstitute these matrices into $\dot{\mathbf{p}}=-\Omega \mathbf{p}+\xi D\mathbf{p}-\Gamma \mathbf{h}$. Projection along and perpendicular to the polar order parameter yieldsAssume $g\ne 0$, $\Gamma >0$, $b>0$ and $p>0$. The flow-alignment angle in planar shear must satisfyFor $\xi =-1$, use the cosine equation: $\theta =\pi /2$ modulo $\pi$, and the displayed tangent is infinite. The radial equation giveswhere the sign distinguishes the angular branches, and only positive right-hand sides are ordered solutions. Linearizing the angular equation gives $\delta \dot{\theta }=-g\xi \sin 2\theta \delta \theta$. For $|\xi |>1$, the stable angular branch therefore hasIts radial relaxation eigenvalue is $-2\Gamma b{p}^2<0$. At $|\xi |=1$ the angular linearization is marginal and $p^2=-a/b$ requires $a<0$. If $g=0$, the shear restriction disappears and any orientation is allowed for $p^2=-a/b>0$. The printed request for dependence only on $a,b,\Gamma ,g$ omits $\xi$: in general the magnitude necessarily depends on it.

Let $h_i=\delta F/\delta p_i$ be the polar molecular field, with the positive functional derivative convention used in the paper. During a pure advective displacement $u_j$, the polar order parameter changes byThe free energy change is $\delta F=\int h_i\delta p_i{d}^{d}r$. Assume periodic boundaries, or boundary conditions that eliminate the surface work, and incompressibility ${\partial }_j{u}_j=0$. Integration by parts in the advective term givesRelabelling the dummy indices in the rotational and flow alignment of a polar order parameter terms givesTherefore $\delta F=\int {\Sigma }_{ij}^p{\partial }_i{u}_j{d}^{d}r$, with the reversible stress of a polar liquid crystalFor a local free-energy density $f(\mathbf{p},\nabla \mathbf{p})$, an explicit choice isIts divergence is exactly $-p_k{\partial }_j{h}_k$, by the Euler-Lagrange equation for $h_k$. A pressure term can be reassigned in an incompressible flow; the chosen representative realizes the printed divergence without that ambiguity. The force density ${\partial }_i{\Sigma }_{ij}^p$ has mechanical power $-\int {\Sigma }_{ij}^p{\partial }_i{v}_j$, the negative of the advective free energy rate, which checks the stress sign.

Use the paper's velocity gradient convention $A_{ij}={\partial }_i{v}_j$, so $\Omega =(A-A^T)/2$ and $D=(A+A^T)/2$. In a rigid motion there is rotation but no rate-of-strain tensor: $D=0$. A polar order parameter attached to the material must rotate exactly with it. In this convention its material derivative is $-\Omega \mathbf{p}$, and its corotational derivative of a polar vector must vanish.

Thus the compensating term is $+\Omega \mathbf{p}$ with coefficient exactly one. This follows from an objective time derivative and is independent of material properties. By contrast, the response to $D\mathbf{p}$ describes deformation and flow alignment of a polar order parameter; its coefficient $\xi$ can depend on the material. The transposed velocity gradient convention is essential to the sign of $\Omega$.

For $q\ne 0$ and $\alpha >0$, each Fourier mode is an Ornstein-Uhlenbeck process with decay rate $r_q=\alpha q^2$. Its explicit Ornstein-Uhlenbeck solution isso the Itô isometry givesTo identify the overdamped Langevin dynamics, choose a friction ${\zeta }_q>0$, a harmonic oscillator potential $V_q={\zeta }_q\alpha q^2{h}_q^2/2$, and force noise $f_q={\zeta }_q{\eta }_q$. The fluctuation-dissipation relation for a Langevin particle then defines $k_B{T}_{\mathrm{eff}}={\zeta }_q{\sigma }^2/2$. Its Boltzmann distribution isThis is the stationary spectrum of a linear fluctuating interface, also following from the equipartition theorem. The numerical prefactor uses exactly the mode-noise normalization given in the paper. The nonzero modes can therefore have a stationary Gaussian distribution even though the unpinned zero mode cannot.

The Neumann boundary conditions give Fourier modes $q=n\pi /L$, including the spatially uniform mode $q=0$. In the frame translating with the mean deposition height, that mode has no restoring force and obeysHere $W_t$ is Brownian motion, independent of the initial height. More generally $⟨h_0(t{)}^2⟩=⟨h_0(0{)}^2⟩+{\sigma }^{2}t$ when these second moments exist. Removing the mean deposition drift does not remove the Gaussian white noise.

The Brownian zero mode of a fluctuating interface has no stationary distribution on the real height axis for ${\sigma }^2>0$. This conclusion does not depend on assuming finite variance: a stationary characteristic function would satisfy $\hat{\mu }(k)=\hat{\mu }(k)e^{-{\sigma }^2{k}^{2}t/2}$, hence vanish for every $k\ne 0$, contradicting its continuity at zero and $\hat{\mu }(0)=1$. Any stationary joint distribution of the whole height would have a stationary zero-mode marginal, which is impossible. Thus the full unpinned height has no Boltzmann equilibrium. Pinning the mean height, or retaining only the nonzero Fourier modes, removes this obstruction. The PDF has ${\dot{h}}_q$; the TeX's ${\hslash }_q$ is a transcription error.

At thermal equilibrium, microscopic reversibility equates the probabilities of a path and its reversed path when both include their Boltzmann distribution initial weights. Denote their endpoint states by $z_1,z_2$, including velocity if needed. Since the Hamiltonian is even under time reversal in classical mechanics,This detailed balance condition and the energy balance for an autonomous Lagrangian yieldThis is the fluctuation-dissipation relation for a Langevin particle: the strength of Gaussian white noise is fixed by the damping and temperature, with $k_B$ the Boltzmann constant.

For an equilibrium coarse-grained variable, the unresolved microscopic states contribute entropy; their statistical weight is encoded in the Helmholtz free energy, rather than in a single microscopic energy. Relative to the same reference measure, ${\rho }_{\mathrm{eq}}(x)\propto e^{-\beta F(x)}$, so microscopic reversibility becomesThis extension assumes an equilibrium coarse-grained description with reversible path statistics; externally driven dynamics need not obey this relation.

For an autonomous Lagrangian, the Euler-Lagrange equation usesApply the chain rule to the Hamiltonian along an arbitrary smooth path, without assuming the unforced Euler-Lagrange equation:Consequently the energy balance for an autonomous Lagrangian isWith explicit time dependence, an additional $-{\partial }_{t}L$ enters $dH/dt$. For ideal Gaussian white noise, the identity is understood through smooth-noise regularization or the Stratonovich chain rule. The original PDF correctly differentiates $L$ with respect to $\dot{x}$ in $H$; the supplied TeX's derivative with respect to $x$, its endpoint $t^2$, and its ordinary derivative of $A$ are transcription errors.

Write $E(t)=\delta A/\delta x(t)$ for the functional derivative of the action. Under time reversal in classical mechanics, $\dot{x}$ changes sign whereas $E$ is even by assumption. Thus the forward and reversed Gaussian white noise histories associated with the same geometric path areThe Onsager--Machlup path probability therefore givesHere the path probabilities are densities conditional on their respective initial states; they are not separately normalized bridges conditioned on both endpoints. The noise-to-path Jacobian determinant must be treated with a consistent discretization: it cancels when it is invariant under time reversal in classical mechanics. This is the usual additive-noise Langevin dynamics convention, including a constant-mass Underdamped Langevin dynamics or the time-reversal invariance of a path Jacobian in midpoint overdamped Langevin dynamics. A completely arbitrary velocity-dependent Lagrangian would require specifying that measure rather than deducing its cancellation solely from the parity of $E$.

The piecewise-linear hat functions have disjoint supports unless their nodes coincide or are neighbours. Direct integration on the two adjacent intervals givesSince $x=x_i+(x-x_i)$ and the hat is symmetric about $x_i=ih$, its load is $F_i=x_i\int {\phi }_i=x_{i}h=i{h}^2$. Thus the reaction-diffusion finite element matrix gives the explicit equationsThe coefficient matrix is a tridiagonal matrix and a positive-definite matrix by the energy identity in part (b), hence nonsingular. A separate spectral verification uses the sine basis with ${\theta }_j=j\pi /(M+1)$ and givesBoth terms are positive for the Dirichlet modes. This also identifies precisely how the diffusion stiffness matrix and reaction mass matrix combine.

The trial space must be a subspace of $H_0^1(0,1)$, not merely of $L^2(0,1)$. The printed assumptions on the basis do not ensure square-integrable weak derivatives, and boundary traces are not defined for arbitrary L2 space functions. For example, $\phi (x)=x^{1/4}(1-x)$ is continuous, vanishes at both endpoints and belongs to $L^2$, but $|{\phi }^{\prime }(x){|}^2\sim x^{-3/2}/16$ is not integrable at zero. Assume the basis belongs to the zero-boundary Sobolev space, as it does for the hat functions in part (c).

Write $u_H={\sum }_{j=1}^M{c}_j{\phi }_j$. Differentiating the restricted energy with respect to $c_i$, or imposing the Galerkin method with test function ${\phi }_i$, givesThe derivative term is the usual diffusion stiffness matrix and the second is the mass matrix. The full matrix is real symmetric, and for every nonzero coefficient vector $c$,because a basis is linearly independent. Thus it is a positive-definite matrix, so the restricted Ritz method has a unique solution and minimum.

Work in the real zero-boundary Sobolev space $X=H_0^1(0,1)$, whose members have square-integrable weak derivatives and zero endpoint trace. Define the bounded bilinear form and linear functionalThe weak formulation is to find $u\in X$ with $a(u,v)=\ell (v)$ for every $v\in X$. Integration by parts recovers the differential equation for a smooth solution; conversely, this identity is the definition of its weak solution.

Use the full $H^1$ norm. The Cauchy-Schwarz inequality gives $|a(u,v)|\le \Vert u{\Vert }_{H^1}\Vert v{\Vert }_{H^1}$ and $|\ell (v)|\le \Vert v{\Vert }_{H^1}/\sqrt{3}$. Also $a(v,v)=\Vert v{\Vert }_{H^1}^2$, so $a$ is a coercive bilinear form with constant one. The Lax-Milgram theorem states that a bounded, coercive bilinear form on a real Hilbert space determines a unique weak solution for every bounded linear functional. It therefore applies here.

The equivalent minimization problem isIf $u$ is the weak solution, symmetry and $a(u,v-u)=\ell (v-u)$ giveThis is positive unless $v=u$, proving a unique global minimizer. Conversely, differentiating $\mathcal{E}(u+tv)$ at $t=0$ gives the weak identity. Thus the unique energy minimizer is exactly the unique weak solution. As a check, the classical representative isIt satisfies both endpoint conditions and $-u^{\prime \prime }+u=x$.

Arrange the interior values into a vector $U$. Represent the five-point Dirichlet Laplacian as a Kronecker sum $L_h=h^{-2}(D\otimes I+I\otimes D)$, where $D$ is the tridiagonal matrix with diagonal $-2$ and adjacent entries one. This Kronecker sum is real symmetric; the real sampled potential has a real diagonal matrix $D_V$. ThusThe generator $-i{H}_h$ is a skew-Hermitian matrix. ConsequentlyThe matrix exponential is a unitary matrix, either by differentiating its product with its adjoint or by unitary diagonalization of a normal matrix. For the two-dimensional discrete L2 norm $\Vert U{\Vert }_h^2=h^2{\sum }_{m,n}|U_{m,n}{|}^2$, this givesApplying the same equation to a difference proves stability of a numerical method with a mesh-independent constant one. This is norm conservation of a semidiscrete Schrödinger equation. It concerns continuous time after spatial discretization; an arbitrary subsequent time integrator need not preserve this stability.

The printed coordinates do not discretize the stated square. For $[-1,1{]}^2$ with $M$ interior points in each direction, use $h=2/(M+1)$, $x_m=-1+mh$, $y_n=-1+nh$, and sample $V$ there. The printed $h=1/(M+1)$ and unshifted coordinates instead describe a grid on $[0,1{]}^2$. The matrix proof is valid for either geometry with its corresponding boundary values, so this transcription-independent statement flaw does not alter the stability conclusion.

For a sufficiently regular solution with zero Dirichlet boundary conditions, differentiate its squared L2 norm. Since $u_t=-i(\Delta u-Vu)$,Integration by parts and the zero boundary trace give $\int \overline{u}\Delta u=-\int |\nabla u{|}^2$. The potential contribution is real because $V$ is real. The quantity in parentheses is consequently real, and its product with $-i$ has zero real part. HenceThis is the energy method for a skew-adjoint evolution. For general L2 space data, the same result follows from the strongly continuous unitary group generated by a self-adjoint realization of $\Delta -V$. A sufficient assumption is a bounded real potential with the Dirichlet Laplacian; arbitrary real potentials require appropriate domain and self-adjointness assumptions. The differentiated proof presumes the regularity needed for the displayed integrals, rather than asserting it for every pointwise real function $V$.

Use $N$ distinct periodic grid points $x_m=m/N$, $0\le m<N$, and $h=1/N$. Values at the two ends of the interval represent the same point, so the endpoint must not be stored twice. With indices interpreted modulo $N$, setThis replaces the Dirichlet matrix by the periodic discrete Laplacian. The normalized Fourier modes $N^{-1/2}{e}^{2\pi ijm/N}$ form an orthonormal basis and have amplification factorsEach satisfies $0<g_j\le 1$ for every $\mu >0$. The Parseval identity therefore gives $\Vert U^{n+1}{\Vert }_h\le \Vert U^n{\Vert }_h$ in the discrete L2 norm, uniformly in the grid and time step. Thus all positive Courant numbers remain stable. The constant mode has $g_0=1$, expressing preservation of the spatial mean; all nonconstant modes decay. The formula with modulo indices also handles small grids, where two neighbours can coincide.

Let $T$ be the $M\times M$ tridiagonal matrix with diagonal two and adjacent entries minus one, incorporating zero Dirichlet boundary conditions. The actual PDF scheme isAn orthonormal basis of eigenvectors of $T$ has components $\sqrt{2/(M+1)}\sin (jm\pi /(M+1))$, and substitution givesThey form an orthonormal basis by the spectral theorem for real symmetric matrices. Hence $I+\mu T$ is invertible for every $\mu >0$, and the amplification eigenvalues are $g_j=(1+\mu {\lambda }_j{)}^{-1}\in (0,1)$. For the discrete L2 norm $\Vert U{\Vert }_h^2=h{\sum }_m|U_m{|}^2$,This bound has constant one independent of time step and mesh, which is the required stability of a numerical method. Therefore every positive Courant number is stable. If a forcing or local-defect sequence is added, the discrete variation-of-constants formula gives a bound by the initial norm plus the sum of those perturbation norms; this also justifies the global error conclusion in part (a).

The original PDF puts both neighbouring values at time $n+1$. Thus the actual scheme is Backward Euler method in time with a centered three-point approximation to the second spatial derivative. The old-time neighbours in the TeX transcription describe a different method.

Let $h=\Delta x$ and $k=\Delta t$. Insert a sufficiently smooth exact solution and divide by $k$. Expanding at $(x_m,t_{n+1})$ givesThe PDE cancels the $u_t-u_{xx}$ term. Therefore the normalized local truncation error is $O(k+h^2)$ andFor compatible smooth data with exact nodal initialization, the stability of a numerical method established below gives global error $O(k+h^2)$ on fixed time intervals, by summing the one-step defects with the contraction bound. Thus no relation between $k$ and $h$ is needed for stability; choosing $k=O(h^2)$ makes both consistency contributions $O(h^2)$.

For algebraic stability of a Runge-Kutta method, the weights must be nonnegative andmust be a positive semidefinite matrix. Here all weights are positive, but direct calculation givesIn particular $e_1^{T}M{e}_1=-1/36<0$; equivalently, its eigenvalues are $0,\pm \sqrt{3}/36$. Therefore the method is not algebraically stable. This does not contradict its A-stability: algebraic stability is a stronger condition designed to guarantee nonlinear B-stability through the Runge-Kutta contractivity identity.

Use the stability function $R(z)=(z^2+6z+12)/(z^2-6z+12)$. Its denominator vanishes only at $3\pm i\sqrt{3}$, both in the open right half-plane. For $z=x+iy$,This is nonnegative whenever $x\le 0$, so $|R(z)|\le 1$ throughout the closed left half-plane. The internal stage equations are also solvable there, since $\det (I-zA)=1-z/2+z^2/12$. Therefore the method is A-stable.

The inequality is strict in the open left half-plane, while $|R(iy)|=1$. Moreover $R(z)\to 1$ as $|z|\to \infty$, so the method is not L-stable. These conclusions follow directly from the definitions of A-stability and L-stability, rather than from the order of the method.

The nodes $0,1/2,1$ have Lagrange interpolation polynomialsThe coefficients obtained from $a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds$ and $b_j={\int }_0^1{\ell }_j(s)ds$ reproduce the given Butcher tableau. In particular, the second row's last entry is $-1/24$, which is missing from the TeX transcription. Thus this is the three-stage Lobatto IIIA method.

One applicable collocation theorem is that $s$-stage Lobatto IIIA collocation has order $2s-2$. Alternatively, the fourth-order conditions for a Runge-Kutta method give a direct verification. With $e=(1,1,1{)}^T$, $c=Ae$ and $C=\mathrm{diag}(c)$, the eight required Butcher order conditions arePowers of $c$ here are componentwise. Substitution satisfies all eight. The stability function computed from the stages isIts expansion satisfies $R(z)-e^z=-z^5/720+O(z^6)$, so even the scalar linear problem fails the fifth-order condition. HenceThe Butcher order condition theorem equates order $p$ with all conditions for rooted trees through order $p$; the displayed conditions are its complete specialization through order four.

Only $-1<\alpha <1$ needs consideration. On the Dahlquist test equation, the amplification polynomial of a multistep method is $P_z(w)=\rho (w)-z\sigma (w)$, where $z=h\lambda$. A-stability requires all amplification roots to have modulus at most one for every $\mathrm{Re}z\le 0$, with simple unit-modulus roots.

Choose the negative real value $z=-2/(1-\alpha )$. The leading coefficient of $P_z$ iswhereas $P_z(1)=-z\sigma (1)=4>0$. Hence $P_z(w)\to -\infty$ as real $w\to +\infty$. The intermediate value theorem supplies a real amplification root $w>1$, so this negative test value is unstable. ThereforeThis is a direct absolute stability obstruction, independent of the Second Dahlquist barrier, which by itself would not rule out the convergent second-order members.