Let $\Delta \rho ={\rho }_p-{\rho }_f>0$. A spherical grain has submerged weightNeglect lift and contact torque, and use a sliding static friction model with normal reaction $R=G$. The inertial quadratic drag iswhere the last equality fixes the convention $\tau ={\rho }_f{U}^2$. Equivalently $U$ is the bed shear velocity and $C_D$ is an effective drag coefficient referred to it. A literal grain-level flow speed can differ from the shear velocity; that conversion must then be absorbed into $C_D$.

Downstream sliding begins when $F_D>{\mu }_{s}R$. Defining the Shields parameter by $\Theta =\tau /(\Delta \rho gd)$, the threshold force balance isThe grain moves downstream above this threshold in the stated sliding model. Real grain motion can instead involve lift, rolling, irregular contacts, or viscous drag; the printed constant belongs to the particular inertial-drag convention and force balance above.

The normal contact force and resisting static friction balance the drag force and submerged weight at impending motion. The right panel uses locally bed-tangent drag force.

For a horizontal source, $m_x=m_0$ throughout. Far downstream $m_z\gg m_0$, so the kinematic plume fluxes obey the vertical pure plume equations to leading order,Set $q=q_{\ast }{f}_0^{1/3}{\hat{s}}^{5/3}$ and $m_z=m_{\ast }{f}_0^{2/3}{\hat{s}}^{4/3}$ with $\hat{s}=s+s_0$. Matching exponents and coefficients gives $(5/3)q_{\ast }=E\sqrt{m_{\ast }}$ and $(4/3)m_{\ast }^2=q_{\ast }$. Consequentlyand the requested literal mass and density-weighted buoyancy fluxes areto leading order. The radius is $b\sim (6\alpha /5)\hat{s}$ and the axial velocity decreases as ${\hat{s}}^{-1/3}$.

The centreline becomes nearly vertical but retains finite horizontal drift:Since $dz/ds=1-O({\hat{s}}^{-8/3})$, there is a finite positive geometric deficit $D={\lim }_{s\to \infty }(s-z)$, and $z=s-D+o(1)$. The equivalent vertical plume virtual origin is therefore located atFor the horizontal point-source forced plume, the virtual source is downstream and below the actual source. The horizontal turbulent round jet first entrains substantial ambient fluid while rising only a little. Thus at the height where it turns upward it already has a finite radius and volume flux; an equivalent pure plume must have begun rising from below to acquire these. The turning displacement, entrained radius divided by its far-field spreading angle, and virtual vertical depth all scale asAt fixed entrainment coefficient this is simply $O(L_J)$. The arclength origin $s=-s_0$ should not be interpreted as a physical height: the curved near field supplies the additional shift $D$. The numerical constants depend on the chosen top-hat plume model and entrainment coefficient; scaling does not fix them.

The inclined forced plume bends from a cubic near-field trajectory toward a vertical asymptote. Dashed far-field spreading lines extrapolate to the equivalent plume virtual origin; the marked location is obtained from the integral model for this illustration, not a universal numerical prediction.

Put $m_0=M_0/{\rho }_0$ and $f_0=F_0/{\rho }_0>0$ under the literal mass convention. Their dimensions are $[m_0]=L^4{T}^{-2}$ and $[f_0]=L^4{T}^{-3}$. Dimensional analysis givesThe jet length compares source momentum with buoyancy. The region near the source behaves as a turbulent round jet; sufficiently far away, the forced plume behaves like a pure plume. With the constant entrainment coefficient retained explicitly, substantial turning occurs over a length of order ${\ell }_J=L_J/\sqrt{E}$, where $E=2\alpha \sqrt{\pi }$. Statements involving $s\ll L_J$ or $s\gg L_J$ normally hold this dimensionless coefficient fixed.

Near the point source, $m\sim m_0$ and the Batchelor entrainment hypothesis gives $q\sim E\sqrt{m_0}s$. ThereforeIntegrating the centreline tangent givesTo leading order the centreline is a straight turbulent round jet at its source inclination. For the horizontal source, the first curvature is the horizontal forced-plume trajectoryFor a vertical source, $\cos {\theta }_0=0$, the trajectory stays vertical rather than developing this cubic transverse displacement. If $F_0=0$, there is no finite buoyancy-induced jet length and the turbulent round jet remains straight.

For a neutrally buoyant horizontal turbulent round jet, a segment's horizontal momentum flux changes only if an external horizontal force or ambient momentum input exists. Neither is present, soPositive buoyancy adds a vertical force, not a horizontal one. Therefore $d{M}_x/ds=0$ remains true for an inclined forced plume even when $F_0\ne 0$; in general $M_x=M_0\cos {\theta }_0$.

The net upward force on a segment is $g({\rho }_0-\rho )\pi b^{2}ds$. Since $F=g({\rho }_0-\rho )\pi b^{2}V$ and $V=|\mathbf{M}|/Q$, vertical momentum balance givesIn kinematic plume fluxes, the equations are $m_x^{\prime }=0$, $m_z^{\prime }=fq/m$, where $m=(m_x^2+m_z^2{)}^{1/2}$. The centreline geometry satisfiesFor $F_0>0$, $Q>0$ beyond the source, and an initially nonvertical turbulent round jet pointing right, $M_z$ increases and the centreline turns toward the vertical. The magnitude of momentum flux is not conserved once the tangent has a vertical component.

It is important to fix the mass flux convention. With literal mass units and the corresponding density-weighted buoyancy flux, the Boussinesq approximation and circular top-hat plume model giveReplacing $\rho$ by ${\rho }_0$ in inertial fluxes is the Boussinesq approximation; the density difference is retained in the buoyancy flux. Define the kinematic plume fluxesIf the symbols $Q,\mathbf{M},F$ instead denote the customary kinematic fluxes, simply omit the factors ${\rho }_0$ in all conversions below. Keeping the two conventions consistent is essential for the jet length and dimensional prefactors.

The Batchelor entrainment hypothesis prescribes inward edge velocityThe constant entrainment coefficient relates ambient inflow to the local axial velocity. A centreline segment of arc length $ds$ entrains volume $2\pi b{u}_{e}ds$. Hence, setting $E=2\alpha \sqrt{\pi }$,The buoyancy flux is conserved because the homogeneous entrained ambient has zero density deficit, and mixing has no buoyancy source or sink. The imposed zero source mass flux with nonzero source momentum flux is a singular point-source idealization, not a finite-radius nozzle with zero emitted fluid. Indeed $V=m/q$ and the plume reduced gravity $f/q$ diverge as $q\to 0$ with nonzero source fluxes. The Boussinesq approximation cannot remain valid arbitrarily close to that point. The top-hat plume model is applied outside the unresolved source region.

Let $\mathcal{V}=H{L}_0$ be the layer volume per unit transverse width. The gravity-current box model imposes $hL=\mathcal{V}$; it approximates the front and scalar balances, rather than an exactly uniform solution of the local momentum equation. For a specified constant front Froude number $\mathrm{Fr}$, the heated particle-laden gravity current satisfiesThe initial data are $L=L_0$, $\varphi ={\varphi }_0$, and $\theta =0$. The gravity-current front condition applies only while $\gamma \varphi -\theta \ge 0$.

For heating without settling, $\varphi ={\varphi }_0$ and $\theta =\beta {\varphi }_{0}t$. Integration of the front equation givesHence the heated gravity-current runout isIt is reached at the neutral-buoyancy time $\gamma /\beta$ within this front model. A finite-momentum current could subsequently coast or lift off: the formula is the maximum predicted by the prescribed gravity-current front condition, which contains no independent front inertia.

For settling without heating, $\theta =0$. Set $K=\mathrm{Fr}\sqrt{g\gamma \mathcal{V}}$, so $\dot{L}=K\sqrt{\varphi }{L}^{-1/2}$. Eliminating time yieldsThus the runout length of a gravity current isThis length is approached as $t\to \infty$, because positive concentration cannot disappear at a finite time under the settling law. If distance means displacement from the removed barrier, the answer is $L_{\max }-L_0$ in either case.

Use the shallow-water approximation, a hydrostatic vertically mixed layer, a deep motionless ambient, and no ambient fluid entrainment. Particle loss changes reduced gravity but changes total layer volume only at the discarded dilute-particle order. Likewise thermal expansion changes the equation of state while the leading Boussinesq approximation retains volume conservation. DefineDepth integration of mass conservation, horizontal momentum balance, particle transport, and the heating law givesHere the integrated excess hydrostatic pressure is $b{h}^2/2$. There is no particle source from the bed. Bed drag and the particle contribution to inertia are neglected at this order. The equivalent variable-buoyancy shallow water equations areThe term $h{b}_x/2$ is necessary when heating or particle loss makes the reduced gravity vary horizontally.

For $v=(h,u,\varphi ,\theta {)}^T$, the quasilinear system $v_t+A(v)v_x=S(v)$ hasIts characteristic polynomial factors asThus the characteristic curves have slopesOn the two repeated contact characteristic curves, the particle and heating laws are $D\varphi =-W_s\varphi /h$ and $D\theta =\beta \varphi$. Writing $c_g=\sqrt{bh}$ and $D_{\pm }={\partial }_t+(u\pm c_g){\partial }_x$, the gravity-wave compatibility relations areFor $h>0$ and $b>0$, gravity-wave right eigenvectors can be taken as $(1,\pm c_g/h,0,0{)}^T$. The contact eigenspace has dimension two: $\delta u=0$ and $b\delta h+(h/2)\delta b=0$, with independent $\delta \varphi ,\delta \theta$. Therefore $A$ has a full set of real eigenvectors.

The system is hyperbolic throughout strict static stability, but not strictly hyperbolic because the contact speed is repeated. At $b=0$ the eigenvalues coalesce and this diagonalizability is lost; when $b<0$ the gravity-wave speeds become complex, reflecting the loss of a statically stable lower layer.

The exact mixture equation of state is $\rho =(1-\varphi ){\rho }_f+\varphi {\rho }_p$. Keeping first-order terms in the small particle volume fraction and thermal expansion givesThe omitted term is ${\rho }_0\varphi \theta$. The downward solid-volume particle deposition flux is $J_v=W_s\varphi$; the corresponding particle mass flux is $J_m={\rho }_p{W}_s\varphi$. Uniform vertical mixing and constant layer depth implyFor $W_s>0$, set $r=W_s/H$. Thenand henceThe heated particle-laden layer loses its stable density stratification when this contrast vanishes. For $\beta >0$, solving for the neutral-buoyancy time givesThe layer is neutral at $T_s$ and statically unstable for $t>T_s$. The subsequent uniformly mixed lower-layer solution cannot represent the resulting overturning.

If $W_s=0$, then $\varphi ={\varphi }_0$, $\theta =\beta {\varphi }_{0}t$, and the continuous limit is $T_s=\gamma /\beta$. If $\beta =0$, the layer stays denser than its surroundings at every finite time and becomes neutral only asymptotically when $W_s>0$; thus $T_s=\infty$. These limiting cases must replace the printed formula when its denominator vanishes.

The central stratified layer supports upward and downward internal gravity waves, whose counterpropagating components form a standing wave. The unstratified regions support evanescent waves, and the upper disturbance decays as $e^{-k(z-2H)}$. Consequently the central layer is a stratified internal-wave guide, rather than a source of propagating energy at infinity.

The largest response occurs at its trapped normal mode frequencies. They can be specified without solving the forced problem. Set the boundary forcing to zero to find a free normal mode. The lower solution is proportional to $\sinh (kz)$, so continuity of fluid pressure and vertical velocity givesIn the central layer writeThe upper Robin boundary condition then gives the exact trapped-mode conditionThe left side increases strictly from zero to infinity, so there is one positive root for each $n$; the frequencies accumulate at zero. This is constructive phase matching after reflection at both ends. The phase shifts from the evanescent waves matter: simply imposing integer half-wavelengths across the stratified layer is generally incorrect.

In the ideal inviscid model, exact resonant forcing has no bounded steady harmonic solution. The undamped normal mode grows secularly under sustained forcing. Weak viscosity or other losses would produce large finite peaks near the displayed frequencies; the linear approximation eventually fails if the disturbance becomes too large.

The illustrated free normal mode has an oscillatory central region and evanescent wave tails. Its lower tail reaches the fixed zero-displacement boundary; its upper tail decays to infinity.

At the change in buoyancy frequency, continuity of vertical velocity and fluid pressure givesThe background mass density is continuous, so there is no interfacial density jump or separate surface restoring force. Since $P=i\omega W^{\prime }/k^2$, these internal-wave transmission across a stratification step conditions are continuity of $W$ and $W^{\prime }$. The buoyancy perturbation need not be continuous, because the background density gradient jumps.

Let $W_0=-i\omega {\eta }_0$ and $\ell =k\sqrt{N_2^2/{\omega }^2-1}$. The lower unstratified layer obeys $W^{\prime \prime }-k^{2}W=0$ and the upper radiation condition givesPropagating these matching data downward givesThe kinematic boundary condition at the moving lower boundary therefore givesThe upper vertical-displacement complex amplitude at $H$ isIn particular, $|{\eta }_H/{\eta }_0|=[{\cosh }^2(kH)+(\ell /k{)}^2{\sinh }^2(kH){]}^{-1/2}$. The unstratified layer attenuates and phase-shifts the transmitted internal gravity wave. The matching height is $H$ in the original PDF; the TeX transcription's $H^{\prime }$ is a typo.

Let $p$ denote the fluid pressure perturbation divided by the reference mass density, and let $b=-g{\rho }^{\prime }/{\rho }_0$ be the buoyancy perturbation. The buoyancy frequency satisfies $N^2=-g{\hat{\rho }}_z/{\rho }_0$. Linearization requires small boundary slope $k|{\eta }_0|\ll 1$, small displacements compared with the vertical scales of the disturbance and background, and advection small compared with the oscillatory acceleration. For a propagating internal gravity wave this includes $|m{\eta }_0|\ll 1$; boundary slope alone is insufficient when $|m|\gg k$. The Boussinesq approximation also requires small relative mass density differences over the region of interest. These conditions must hold for the resulting disturbance, including any amplification by resonance.

The inviscid Linearized Boussinesq equations and linear kinematic boundary condition areWrite the velocity field as $(u,w)=\mathrm{Re}\{(U(z),W(z))e^{i(kx-\omega t)}\}$. Incompressibility and horizontal momentum balance giveThe vertical momentum balance therefore yieldsFor a real vertical wavenumber $m$, this gives the dispersion relationWhen $0<\omega <N$, define $\ell =k\sqrt{N^2/{\omega }^2-1}$. The radiation condition selects $m=-\ell$, because energy must leave the boundary upward. The boundary-forced internal gravity wave isIts wavevector is $\mathbf{K}=(k,m)$, its phase velocity is ${\mathbf{c}}_p=\omega \mathbf{K}/(k^2+m^2)$, and its group velocity isThus the phase velocity points down and right, while the group velocity points up and right. They are perpendicular. Internal-wave polarization makes the particle motion an oscillation along the group velocity direction, with zero first-order mean transport. The constant-phase lines of an internal gravity wave are also parallel to the group velocity. Their inclination $\vartheta$ above the horizontal satisfies $\sin \vartheta =\omega /N$.

When $\omega >N$, set $\kappa =k\sqrt{1-N^2/{\omega }^2}$. Boundedness at infinity selects the evanescent waveThe horizontal and vertical velocity components are in quadrature: fluid particles describe small ellipses, and the response decays over ${\kappa }^{-1}$. There is no upward time-averaged energy flux, because $P$ and $W$ are in quadrature. A real vertical group velocity is not defined for this evanescent wave. The pattern travels horizontally with phase velocity $\omega /k$.

For $N=0$ the response is the evanescent wave with $\kappa =k$; the ratios $\omega /N$ should not be used. At the cutoff $\omega =N>0$, the bounded harmonic solution has $W=-i\omega {\eta }_0$, $U=0$: it neither decays nor has nonzero upward group velocity. It is the limiting cutoff response, rather than a localized radiating disturbance. A uniform nonzero buoyancy frequency in an infinitely deep Boussinesq approximation is itself a local idealization of the background mass density.

The arrows distinguish the phase velocity, group velocity, and oscillatory particle motion of the internal gravity wave. The right panel shows the decay envelope and particle ellipses of the evanescent wave.

Within this analytic quartic Landau theory, equilibrium of a nonconserved scalar order parameter requires $f^{\prime }=0$. At a continuous phase transition, its restoring curvature vanishes, so $f^{\prime \prime }=0$. A nonzero $f^{\prime \prime \prime }$ would leave a cubic leading term in the Taylor expansion about that state and prevent it from being a local minimum. Thus a stable critical endpoint of a quartic Landau free energy requires all three derivatives to vanish, with $f^{\prime \prime \prime \prime }>0$. The condition $f^{\prime \prime }=0$ alone can instead identify a spinodal point.

For $f=\bar{b}{\lambda }^4+\bar{c}{\lambda }^3+\bar{a}{\lambda }^2+\bar{d}\lambda$, $\bar{b}>0$, the successive equations areThey giveAt these values, the free energy is exactly $f({\lambda }_c)+\bar{b}(\lambda -{\lambda }_c{)}^4$, proving stability. With no electric field, $\bar{d}=0$, so necessarily $\bar{c}=0$, followed by ${\lambda }_c={\bar{a}}_c=0$.

For $\bar{c}<0$ and $\chi >0$, the electric-field coupling to nematic order supplies the required negative ${\bar{d}}_c$ atTuning the quadratic coefficient to ${\bar{a}}_c$ simultaneously reaches the field-induced critical endpoint of the isotropic-nematic transition. The PDF's $\tilde{c}$ in the last clause is undefined; it is interpreted here as the previously defined $\bar{c}$. The endpoint is a prediction of the stable quartic, small-field mean-field approximation; higher-order terms or a field beyond the expansion's validity can shift it.

A scalar coupling linear in the nematic order parameter requires a second-rank anisotropy. A single electric field supplies no such tensor, while $E_i{E}_j$ does. Its isotropic part gives $E^2\mathrm{Tr}Q=0$. Hence the lowest joint order in $E$ and $Q$ of the electric-field coupling to nematic order isTerms depending only on $E$ do not affect minimization over $Q$. With $\vartheta$ the angle between the nematic director and the electric field,For $\chi \lambda >0$, the lowest free energy occurs at ${\cos }^2\vartheta =1$, so the nematic director aligns along either $+\mathbf{E}$ or $-\mathbf{E}$. Then $E_i{Q}_{ij}{E}_j=\lambda E^2$, givingUnlike the linear composition term in a conserved binary fluid mixture, this term changes the equilibrium nematic order parameter: its spatial integral is not fixed.

For the stated uniaxial nematic order normalization, the nematic director is an eigenvector of $Q$ with eigenvalue $\lambda$, and the two perpendicular eigenvalues are $-\lambda /2$. ThereforeUsing precisely the Landau-de Gennes free energy convention in the preceding solution givesIf the quartic invariant were instead normalized as $b\mathrm{Tr}{Q}^4/4$, its coefficient would be $\bar{b}=9b/32$; the physical predictions are unchanged after redefining $b$.

For a nonzero global minimizer of the free energy, compare opposite values of the scalar nematic order parameter:The lower one has $\bar{c}\lambda <0$, soThe claim concerns the stable ordered phase, not every metastable stationary point. For $\bar{c}=0$ the two signs are degenerate. For $\bar{c}\ne 0$, equality of the ordered and isotropic free energies, together with stationarity, gives ${\lambda }_t=-\bar{c}/(2\bar{b})$ and ${\bar{a}}_t={\bar{c}}^2/(4\bar{b})>0$. The finite jump ${\lambda }_t$ exhibits the first-order phase transition driven by the cubic invariant.

The nematic order parameter is a symmetric second-rank tensor and a traceless second-rank tensor describing orientational anisotropy. For molecular unit axes $\mathbf{s}$, $Q_{ij}$ is proportional to $⟨s_i{s}_j-{\delta }_{ij}/3⟩$. It vanishes in an isotropic phase and is unchanged by head-tail reversal $\mathbf{s}\mapsto -\mathbf{s}$, unlike a polar order parameter. For uniaxial nematic order, the nematic director is thus an unoriented axis.

One convention for the Landau-de Gennes free energy through fourth order, respecting rotational symmetry, isFor a traceless second-rank tensor in three dimensions, the Cayley-Hamilton theorem gives $\mathrm{Tr}{Q}^4=(\mathrm{Tr}{Q}^2{)}^2/2$, so there is only one independent quartic rotational invariant of a symmetric traceless tensor. The only possible linear scalar is $\mathrm{Tr}Q=0$; without an external anisotropy no linear term survives.

The cubic rotational invariant of a symmetric traceless tensor distinguishes prolate and oblate forms of uniaxial nematic order. Head-tail reversal leaves $Q$ unchanged and does not impose $Q\mapsto -Q$. Thus a cubic term is allowed in three dimensions and generically produces a first-order phase transition, rather than a symmetry-enforced continuous onset. In two dimensions a symmetric traceless second-rank tensor has eigenvalues $s,-s$, hence $\mathrm{Tr}{Q}^3=0$: the cubic invariant vanishes identically in two dimensions. This is a statement about the Landau free energy; fluctuations in two dimensions require a separate treatment.

Take $b>0$ for a stable quartic Landau free energy. The cubic-term removal in a conserved quartic Landau free energy uses $\varphi =\psi -c/(3b)$ and gives, apart from a constant,For a closed binary fluid mixture, the compositional order parameter has fixed integral. Adding $d\int \varphi$ therefore adds the same constant to every allowed state. Equivalently, in the common-tangent construction for phase coexistence, a linear term changes the tangent slope, hence the chemical potential, by $d$, but leaves the coexistence compositions unchanged. This is why linear composition bias leaves coexistence compositions unchanged; it would not mean an unchanged composition at a fixed externally imposed chemical potential.

In the shifted variable, the binodal compositions are ${\psi }_{\pm }=\pm \sqrt{-a^{\prime }/b}$ for $a^{\prime }<0$. They merge atThe continuous phase transition therefore has the same scalar quartic Landau theory as the symmetric binary fluid mixture, with a shifted critical composition and control parameter. In the mean-field approximation, the composition difference vanishes with order-parameter critical exponent $1/2$. This argument presumes the critical composition lies in the physically accessible composition range; conservation affects the dynamics, not this static phase coexistence construction.

At $\xi =0$, the polar order parameter obeysFor $a<0$, $b>0$, the magnitude relaxes to $\sqrt{-a/b}$ while the direction rotates continuously. Thus the likely ordered behaviour is tumbling rather than a fixed alignment.

More generally, for $|\xi |<1$ and $g\ne 0$, $\xi \cos 2\theta -1$ never vanishes, so the tumbling of a polar order parameter persists. The polar order parameter completes a full $2\pi$ rotation inobtained by integrating $dt=2|d\theta |/(|g|(1-\xi \cos 2\theta ))$ over a full rotation. An unoriented nematic director repeats after half this period, whereas a polar order parameter distinguishes opposite directions. For $\xi \ne 0$ its magnitude can also oscillate. If the polar order parameter decays to zero, its orientation is no longer a physical observable; the tumbling conclusion concerns a persistent ordered state.

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.