Pith. sign in

REVIEW 3 major objections 3 minor 42 references

Discontinuous phase transition in chemotactic aggregation with density-dependent pressure

T0 review · 3 major / 3 minor · reviewed 2026-08-14 · deepseek-v4-flash

Pith's one-line read The aggregation transition in a two-dimensional chemotaxis model with density-dependent pressure is discontinuous: the order parameter jumps, the Lyapunov landscape has two coexisting minima, and direct simulation exhibits hysteresis.

desk verdict Careful Lyapunov-functional analysis of a chemotaxis model claims a discontinuous aggregation transition with hysteresis, but the numerical evidence rests on a single low-order mode truncation that needs a convergence check. read the letter →

arxiv 1908.08842 v2 pith:4FKEELJP submitted 2019-08-23 cond-mat.stat-mech

classification cond-mat.stat-mech MSC 35Q9292C1735B36 PACS 64.60.-i
keywords chemotaxisPatlak-Keller-Segeldensity-dependentpressureLyapunovfunctionaldiscontinuousphasetransitionhysteresisFourier-Besselexpansionaggregation
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper studies a two-dimensional Patlak-Keller-Segel model in which chemotactic attraction is balanced by a density-dependent pressure that is quadratic in the local density. Its central claim is that the transition from a uniform population to an aggregated one is discontinuous: the order parameter $A = \max\rho - \min\rho$ jumps at a threshold, the Lyapunov landscape develops two coexisting minima, and direct numerical simulation shows hysteresis, so an existing aggregate survives when the attraction strength drops below the linear instability threshold $f_0^{\mathrm{unstable}}\approx2.195$. A two-mode Fourier-Bessel approximation shows why: the quadratic pressure term couples modes, and once the lowest mode is excited, the coupling can create a metastable aggregate branch before the homogeneous solution loses stability. If the claim is right, bacterial aggregation in confined geometries is an abrupt, history-dependent event rather than a smooth condensation.

What carries the argument

The central object is the Lyapunov functional $W$ of the system, equation (3), which decreases along every trajectory and is interpreted as a free energy. The argument is carried by expanding $\rho$ and $c$ in Fourier-Bessel modes $J_p(j'_{pm} r/l)\cos[p(\theta-\eta_{pm})]$ on the disk, where $j'_{pm}$ are zeros of the derivative of the Bessel function; this makes $W_{\mathrm{st}}$ a sum over modes plus a coupling term. The crucial identity is the mode-coupling coefficient $I_{112}=\int_0^1 x J_1^2(j'_{11}x)J_2(j'_{21}x)\,dx\approx0.0474258$, which appears in the cubic term $\kappa I_{112}z_{11}^2 z_{21} R_{11}^2 R_{21}\cos 2(\eta_{11}-\eta_{21})$ in Eq. (8). This term is what produces a second local minimum before the homogeneous solution becomes linearly unstable, and the same Bessel orthogonality relations are used to derive the quadratic form for $\kappa=0$ and the radial-symmetry quartic analysis.

What would settle it

Run the full PDE Eq. (1) with many Fourier-Bessel modes or a high-resolution finite-element scheme across $f_0=2.15$ through $2.26$ and measure the order parameter $A$ on slow up- and down-ramps; if the aggregate dissolves immediately below $f_0^{\mathrm{unstable}}\approx2.195$, or if $A$ grows continuously instead of jumping, the discontinuous-transition claim fails.

Watch

Extended reading notes

Core claim

On the paper's own terms, the main discovery is that adding a quadratic term to the effective pressure changes the character of the aggregation instability expected from the classical PKS model. For $\kappa=D_1/D_0=0$, each Fourier-Bessel mode contributes independently to the Lyapunov functional and the first excited mode $j'_{11}$ produces a jump with no hysteresis. For $\kappa>0$, the cubic coupling term $\propto R_{11}^2 R_{21}$ makes the landscape bistable: the aggregate branch coexists with the uniform branch, so minimization of $W_{\mathrm{st}}$ gives a crossing of branches and a discontinuous jump of $A$, while spectral simulation of the PDE shows that an aggregate formed at $f_0=2.26$ dissolves only at $f_0\lesssim2.15$, well below $f_0^{\mathrm{unstable}}\approx2.195$. The authors also show the boundary matters: an aggregate forms near the disk edge, and under imposed radial symmetry the same model gives a continuous transition with order parameter growing linearly with slope $\propto \kappa^{-1}$.

Load-bearing premise

The argument depends on truncating the infinite Fourier-Bessel expansion to a small number of modes ($P=M=3$ in the numerics, two modes in the analytic picture), and the paper never checks that the jump and hysteresis survive when more modes are included.

Editorial extensions

If this is right

  • Aggregation is not gradual: as chemotactic strength $f_0$ rises through the threshold, a compact aggregate forms at the boundary with a sudden jump in the order parameter $A$.
  • The transition is hysteretic: an aggregate created at high $f_0$ remains intact down to about $f_0\simeq2.15$, so reversing the control parameter does not retrace the forward path.
  • The quadratic pressure term is essential: for linear pressure ($\kappa=0$) the mode analysis gives a coexistence-free jump, while $\kappa>0$ creates bistability and hysteresis.
  • If aggregation is forced to be radial about the disk center, the transition becomes continuous with $A\propto \kappa^{-1}$ above threshold, so boundary geometry can change the order of the transition.
  • The Lyapunov functional can serve as a free energy; a derivative of the minimized functional with respect to disk area acts like a thermodynamic pressure, suggesting a thermodynamic reading of the aggregation transition.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the discontinuity survives the full infinite-mode limit, then a microfluidic or swarm experiment that slowly ramps chemoattractant up and down should observe a hysteresis loop in population clumpiness, not just in the model; this is directly testable with bacterial or synthetic chemotactic swimmers.
  • The paper's own lack of convergence checks for the truncation sizes $P$ and $M$ means the strongest plausible threat is that high-order modes round off the jump or destabilize the aggregate branch; a high-resolution spectral or finite-element study of Eq. (1) is the natural check.
  • For active-matter thermodynamics, the two-branch structure of $W_{\min}$ suggests defining a latent-like quantity from the derivative of the minimized functional, which could support a Maxwell equal-area construction for the aggregation transition.
  • Because Appendix A derives the Lyapunov functional for an arbitrary $\rho$-dependent pressure, the same two-mode criterion could classify supercritical versus subcritical aggregation for other pressure virial forms, not just the quadratic truncation studied here.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 3 minor

Summary. The paper studies a two-dimensional Patlak-Keller-Segel model with density-dependent pressure, where the effective pressure is truncated at quadratic order in density. The authors construct an exact Lyapunov functional and analyze it by expanding the density and chemoattractant fields in Fourier-Bessel modes on a disk. For the linear-pressure case (κ=0) the functional decouples into independent modes and the transition to aggregation is discontinuous but without hysteresis. For nonlinear pressure (κ>0), a two-mode approximation reveals bistability, and numerical minimization of the truncated Lyapunov functional together with direct spectral simulations, both using P=M=3 modes, show a discontinuous transition with hysteresis. A radially symmetric restriction (p=0) is also analyzed and found to yield a continuous transition. The paper concludes that the full aggregation transition is discontinuous and interprets the Lyapunov functional as a free energy.

Significance. If the central claim holds, the paper provides an interesting example of a discontinuous, hysteretic aggregation transition in a chemotaxis model with a mechanical regularization, and it demonstrates the usefulness of Lyapunov-functionals as free-energy analogues for active matter. The derivations in the appendices are careful and transparent: the Lyapunov functional is derived from the equations rather than assumed, and the linear stability threshold is a derived quantity. The numerical methods are clearly described and the hysteresis signal is visually consistent. However, the significance is tempered by the fact that the main quantitative evidence for the full-PDE claim rests entirely on a low-order Fourier-Bessel truncation (P=M=3) without any convergence test, and the paper's own radial-symmetry analysis shows that the transition order can change when the allowed mode structure is changed.

major comments (3)
  1. [§III, Figs. 2 and 3] The central claim that the aggregation transition is discontinuous in the full two-dimensional PDE is supported only by calculations with P=M=3 modes. No convergence study in P and M is reported, and the two numerical methods (gradient-descent minimization of the Lyapunov functional and the spectral simulation) are not independent with respect to truncation because they project onto the same low-order subspace. Since the aggregate forms near the boundary, where the first unstable mode J_1(j'_{11}r/l) peaks, an adequate representation may require many higher-order radial and angular modes. Please add a systematic convergence test (e.g., increasing P and M to 4, 5, or 6) for both the Lyapunov-minimization branches and the hysteresis loop, or provide a rigorous error bound showing that neglected modes do not alter the stability of the clustered state or the shape of the Lyapunov landscape.
  2. [§IV, Eq. (13)] The paper's own radially symmetric analysis shows that when the allowed mode structure is restricted to p=0 and m=1,2, the transition becomes continuous, with the order parameter increasing linearly near threshold. This is explicitly acknowledged in the Discussion as a case where 'the transition behavior may change if we suppress aggregation at the boundary.' This example demonstrates that the order of the transition is sensitive to the modal truncation, which reinforces the need for a convergence check in the full p>0 case. Absent such a check, the discontinuous-transition claim is not established beyond the truncated model.
  3. [Appendix C, Eq. (C3)] The dispersion relation is misprinted: the line reads 'η^2 + [ ... ] + k^2(...)(...)' with no linear term in η and no '= 0'. This makes it difficult to verify the derivation of the stability condition in Eq. (C4). Please correct the equation.
minor comments (3)
  1. [Fig. 2 caption] The caption states the minimization uses P=M=3, but the rendered panel (c) appears to include a legend with 'P=M=5' alongside '3'. If panel (c) shows mode coefficients for two different truncations, the text and caption should state this explicitly; if not, the stray 'P=M=5' label should be removed.
  2. [Appendix C, Eq. (C2)] The matrix notation in the linearized equation is slightly unconventional; the diffusion term is written as a matrix acting on ∇^2 of the perturbation vector. This is understandable, but a sentence explaining the matrix structure would improve clarity.
  3. [Sec. II, Eq. (8)] The two-mode approximation is introduced as an illustration, and the conclusion is explicitly limited to that truncation. This is fine, but the transition from 'the answer is negative, at least within the two-mode approximation' to the later claim about the full system would benefit from a clearer statement that the full-system conclusion relies on the numerical results.

Circularity Check

0 steps flagged · score 2.0 of 10

No significant circularity: the Lyapunov functional, stability threshold, and hysteresis are derived from the model equations rather than fitted or assumed; the sole self-citation sets a parameter value (χ0=4) and is not load-bearing.

full rationale

The paper's central claim is that the aggregation transition in the PKS model with quadratic density-dependent pressure is discontinuous. Each load-bearing ingredient is derived rather than assumed. The Lyapunov functional W in Eq. (3) is constructed in Appendix A from the governing equations (1), with dW/dt ≤ 0 proved from the dynamics and boundary conditions, so minimizing W is not an input equivalent to the target claim. The linear instability threshold f0_unstable ≈ 2.195 follows from Eq. (9), which is derived from the coefficient of R11^2 in the mode-expanded functional and is explicitly checked against the linear stability analysis in Appendix C; it is not a fitted parameter. The two-mode expression Eq. (8) is a truncation of the same derived functional, and the multi-mode gradient-descent minimization in Sec. III A and the spectral simulation of Eq. (1) in Sec. III B are independent numerical routes to the hysteresis and bistability shown in Figs. 2 and 3. The only self-citation, Ref. [27], is used to set χ0=4; this is a parameter choice that rescales the threshold but does not force the discontinuous transition. The shared Fourier-Bessel truncation P=M=3 is a convergence risk (the paper does not test higher modes, and the radial-symmetry calculation in Sec. IV shows sensitivity to allowed modes), but that is a correctness or robustness concern, not circularity: no fitted quantity is renamed as a prediction and no result is assumed by construction. Accordingly, no circular step is exhibited; the score of 2 reflects only the minor, non-load-bearing self-citation for the parameter value χ0=4.

Assumptions & free parameters 0 free parameters · 3 assumptions · 0 invented entities

No free parameters are fitted to data; the model parameters are set to fixed values following Ref. 27 or varied as a control parameter. The main load-bearing assumptions are the quadratic pressure truncation, the boundary conditions, and the low-rank mode truncation used in both the analysis and the numerics.

assumptions (3)
  • domain assumption The effective pressure is truncated at quadratic order, Π(ρ) = D0 ρ + D1 ρ²/2.
    This defines the model; the paper states this truncation in Eq. (2) and uses it throughout.
  • domain assumption The disk boundary conditions are zero-flux Neumann conditions with total mass conservation.
    Stated in Sec. II. This is a modeling choice for a Petri-dish setup.
  • ad hoc to paper Truncation at P=M=3 captures the order of the phase transition.
    The numerical analysis in Sec. III uses only 3 modes; no convergence test is presented, so the validity of this assumption determines whether the discontinuous transition persists in the full model.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Discontinuous phase transition in chemotactic aggregation with density-dependent pressure." pith.science (2026). https://pith.science/paper/4FKEELJP

@misc{pith2026190808842,
  author       = {Pith},
  title        = {Pith review of: Discontinuous phase transition in chemotactic aggregation with density-dependent pressure},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/4FKEELJP}},
  note         = {Machine review of arXiv:1908.08842}
}
read the original abstract

Many small organisms such as bacteria can attract each other by depositing chemical attractants. At the same time, they exert repulsive force on each other when crowded, which can be modeled by effective pressure as an increasing function of the organisms' density. As the chemical attraction becomes strong compared to the effective pressure, the system will undergo a phase transition from homogeneous distribution to aggregation. In this work, we describe the interplay of organisms and chemicals on a two-dimensional disk with a set of partial differential equations of the Patlak-Keller-Segel type. By analyzing its Lyapunov functional, we show that the aggregation transition occurs discontinuously, forming an aggregate near the boundary of the disk. The result can be interpreted within a thermodynamic framework by identifying the Lyapunov functional with free energy.

Figures

Figures reproduced from arXiv: 1908.08842 by the authors.

Figure 2
Figure 2. FIG. 2. (a) Aggregational part of the Lyapunov functional, [PITH_FULL_IMAGE:figures/full_fig_p003_2.png] view at source ↗
Figure 3
Figure 3. FIG. 3. Direct simulation of Eq. (1) using the spectral [PITH_FULL_IMAGE:figures/full_fig_p004_3.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

42 extracted references · 36 canonical work pages

  1. [1]

    (C3) For η to be negative, the following inequality must be met: k2 > ( 1 + D1 D0 ρ0 ) −1 f0χ0 ν0D0 − g0 ν0 , (C4) which is the stability condition

    + g0 + k2ν0 ] + k2(D0ρ0 + D1ρ2 0)(g0 + k2ν0). (C3) For η to be negative, the following inequality must be met: k2 > ( 1 + D1 D0 ρ0 ) −1 f0χ0 ν0D0 − g0 ν0 , (C4) which is the stability condition. 10 Appendix D: Two-mode approximation under radial symmetry When we take ρ ≈ ρ0 + Q01J0(j′ 01r/l) + Q02J0(j′ 02r/l), the corresponding Lyapunov functional becomes...

  2. [2]

    (D1) does not contain R02 and R01R02 due to the orthogonality between normal modes

    (D3) We note that Eq. (D1) does not contain R02 and R01R02 due to the orthogonality between normal modes. Therefore, when R02 ∼ O ( R01 2) around the minimum, the lowest contribution from R02 is of O ( R4 01 )

  3. [3]

    J. A. Shapiro, Sci. Am. 258, 82 (1988)

  4. [4]

    Ma and J

    M. Ma and J. W. Eaton, Proc. Natl. Acad. Sci. U.S.A. 89, 7924 (1992)

  5. [5]

    J. A. Shapiro, BioEssays 17, 597 (1995)

  6. [6]

    B. J. Crespi, Trends Ecol. Evol. 16, 178 (2001)

  7. [7]

    Detrain and J.-L

    C. Detrain and J.-L. Deneubourg, Phys. Life Rev. 3, 162 (2006)

  8. [8]

    P. S. Stewart and J. W. Costerton, Lancet 358, 135 (2001)

Show all 42 references
  1. [9]

    P. R. Secor, L. A. Michaels, A. Ratjen, L. K. Jennings, and P. K. Singh, Proc. Natl. Acad. Sci. U.S.A. 115, 10780 (2018)

  2. [10]

    C. S. Patlak, Bull. Math. Biophys. 15, 311 (1953)

  3. [11]

    E. F. Keller and L. A. Segel, J. Theor. Biol. 26, 399 (1970)

  4. [12]

    Nanjundiah, J

    V. Nanjundiah, J. Theor. Biol. 42, 63 (1973)

  5. [13]

    Childress and J

    S. Childress and J. K. Percus, Math. Biosci. 56, 217 (1981)

  6. [14]

    Stevens, SIAM J

    A. Stevens, SIAM J. Appl. Math. 61, 183 (2000)

  7. [15]

    Hillen and K

    T. Hillen and K. J. Painter, J. Math. Biol. 58, 183 (2009)

  8. [16]

    Gamba, D

    A. Gamba, D. Ambrosi, A. Coniglio, A. de Candia, S. Di Talia, E. Giraudo, G. Serini, L. Preziosi, and F. Bussolino, Phys. Rev. Lett. 90, 118101 (2003)

  9. [17]

    Kowalczyk, J

    R. Kowalczyk, J. Math. Anal. Appl. 305, 566 (2005)

  10. [18]

    Kowalczyk and Z

    R. Kowalczyk and Z. Szyma´ nska, J. Math. Anal. Appl. 343, 379 (2008)

  11. [19]

    Sinhuber and N

    M. Sinhuber and N. T. Ouellette, Phys. Rev. Lett. 119, 178003 (2017)

  12. [20]

    Wittkowski, A

    R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo, and M. E. Cates, Nat. Commun. 5, 4351 (2014)

  13. [21]

    Tiribocchi, R

    A. Tiribocchi, R. Wittkowski, D. Marenduzzo, and M. E. Cates, Phys. Rev. Lett. 115, 188302 (2015)

  14. [22]

    S. C. Takatori, W. Yan, and J. F. Brady, Phys. Rev. Lett. 113, 028103 (2014)

  15. [23]

    S. C. Takatori and J. F. Brady, Phys. Rev. E 91, 032117 (2015)

  16. [24]

    Horstmann, Colloq

    D. Horstmann, Colloq. Math. 87, 113 (2001)

  17. [25]

    Biler, Adv

    P. Biler, Adv. Math. Sci. Appl. 8, 715 (1998)

  18. [26]

    Calvez and L

    V. Calvez and L. Corrias, Commun. Math. Sci. 6, 417 (2008)

  19. [27]

    Fatkullin, Nonlinearity 26, 81 (2013)

    I. Fatkullin, Nonlinearity 26, 81 (2013)

  20. [28]

    K. G. Petrosyan and C.-K. Hu, Phys. Rev. E 89, 042132 (2014)

  21. [29]

    S. K. Baek and B. J. Kim, Sci. Rep. 7, 8909 (2017)

  22. [30]

    M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993)

  23. [31]

    Fanelli, C

    D. Fanelli, C. Cianci, and F. Di Patti, Eur. Phys. J. B 86, 142 (2013)

  24. [32]

    Madzvamuse, H

    A. Madzvamuse, H. S. Ndakwo, and R. Barreira, J. Math. Biol. 70, 709 (2015)

  25. [33]

    Gambino, M

    G. Gambino, M. Lombardo, and M. Sammartino, Acta Appl. Math. 132, 283 (2014)

  26. [34]

    M. L. Boas, Mathematical Methods in the Physical Sci- ences, 3rd ed. (Wiley, Hoboken, NJ, 2006)

  27. [35]

    R. M. Ford, B. R. Phillips, J. A. Quinn, and D. A. Lauf- fenburger, Biotechnol. Bioeng. 37, 647 (1991)

  28. [36]

    Saragosti, V

    J. Saragosti, V. Calvez, N. Bournaveas, B. Perthame, A. Buguin, and P. Silberzan, Proc. Natl. Acad. Sci. U.S.A. 108, 16235 (2011)

  29. [37]

    M. E. J. Newman, Computational Physics (CreateSpace Independent, San Bernardino, CA, 2013). 11

  30. [38]

    Sittner, E

    A. Sittner, E. Bar-David, I. Glinert, A. Ben-Shmuel, S. Weiss, J. Schlomovitz, D. Kobiler, and H. Levy, PloS one 12, e0186613 (2017)

  31. [39]

    Bonazzi, V

    D. Bonazzi, V. L. Schiavo, S. Machata, I. Djafer- Cherif, P. Nivoit, V. Manriquez, H. Tanimoto, J. Husson, N. Henry, H. Chat´ e,et al. , Cell 174, 143 (2018)

  32. [40]

    C. M. Bender, B. K. Berntson, D. Parker, and E. Samuel, Am. J. Phys. 81, 173 (2013)

  33. [41]

    A. J. T. M. Mathijssen, F. Guzm´ an-Lastra, A. Kaiser, and H. L¨ owen, Phys. Rev. Lett. 121, 248101 (2018)

  34. [42]

    A. J. T. M. Mathijssen, R. Jeanneret, and M. Polin, Phys. Rev. Fluids 3, 033103 (2018)

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.