Papers

We summarize recent developments in the theory of the Benjamin-Ono equation, a well-known long-wave asymptotic model for internal waves in deep water, focusing on results obtained by means of an explicit formula for the solution of the Cauchy problem in terms of given initial data. We give a full description of this formula, and provide simple examples, with the aim of bringing it into the realm of applied mathematics. We also show how the formula simplifies further upon restriction to rational initial data. We then show how the formula and its rational restriction explain in great detail, and in a strikingly simple fashion, the asymptotic behavior of solutions of the Benjamin-Ono equation in the small-dispersion and long-time limits. In the setting of long-time asymptotics, we prove some new results related to soliton resolution yielding pointwise convergence with a concrete decay rate under the assumption of rational initial data. The first half of the paper is a stand-alone general survey and user's guide, while the appendices constituting the second half are intended for readers who want to explore the topics in greater depth, and it includes the proofs of our new results.

We review the discrete theory of differential-geometric Poisson brackets as introduced by B. A. Dubrovin, and we present a complete proof of their characterisation in the non-degenerate case. We use this characterisation to show several novel examples of such structures, providing a complete classification for non-degenerate four-dimensional differential-geometric Poisson brackets under the additional assumption of the Lie algebra to be of quasi-Frobenius type. We also present a complete classification of quasi-Frobenius filiform N-graded Lie algebras, including two families in arbitrary dimensions.

It is widely expected that the von Neumann entanglement of a local operator grows at most logarithmically in integrable many-body systems in infinite volume. We give a counterexample in a four-state brickwork circuit whose gate is a permutation matrix solving the constant Yang–Baxter equation. The von Neumann operator entropy grows as $(\log2)\sqrt{t/\pi}+O(\log t)$, whereas fixed-index Rényi entropies grow linearly below index one and logarithmically above it. At fixed relative Hilbert–Schmidt error below one, matrix product operator simulations require a bond dimension growing at least as $\exp(c\sqrt t)$ for some $c>0$.

The Kaleidoscope Yang-Baxter equation of Qiu, Guan, and Yu (2026) is the consistency condition of multiple scattering in Gaudin's kaleidoscope models. At an even order $N$ it involves two matrices: a shift, and a square-zero matrix with one complex parameter. They conjectured that every chained word, the square-zero matrix alternating with integer powers of the shift, vanishes once the number of powers reaches $N/2$, verified through order ten. We prove it for every even order and every parameter where the matrices exist. The bound is sharp. One step below it, a word is nonzero exactly when no factor's power times its position is divisible by half the order, and then has rank two. The proof rests on one flag that every factor lowers by a step. The algebra the two matrices generate is a monotone-path algebra: its basis paths multiply only while their level keeps direction. It is independent of the parameter and has dimension $N^2-N+2$. Its radical, generated by the square-zero matrix, is nilpotent of index one more than half the order, and its representation type is infinite. Four conditions over an arbitrary field reproduce the threshold and the radical. None can be dropped, and the root-of-unity model is one instance.

We investigate the recently introduced polygonal hypergeometric functions [arXiv:2502.12127, arXiv:2507.01904], which are conjectured to evaluate multipoint parametric one-loop conformal integrals in arbitrary dimensions. These functions can be systematically constructed using a diagrammatic algorithm that operates in terms of simple planar shapes. Given an $n$-point polygonal power series, we formulate the corresponding GKZ hypergeometric system, which admits an explicit diagrammatic interpretation. In particular, the transianic matrix, which encodes all the relevant diagrammatic data, determines the corresponding toric matrix along with a basis of the lattice of relations, thereby singling out a finite subsystem of $n(n-3)/2$ toric equations. This finite subsystem consists of second- and third-order equations, whereas the full toric family also includes higher-order ones. We illustrate the general procedure with several prominent examples of polygonal functions, including the fourth Appell and the Srivastava–Daoust functions.

This work aims to investigate the long-time asymptotic behavior of solutions to the Cauchy problem for the Heisenberg ferromagnet equation with the Schwartz initial data. Utilizing the gauge transformations, spectral analysis and the inverse scattering method, we prove that the solutions to the Heisenberg ferromagnet equation can be expressed in terms of the solutions to a matrix Riemann-Hilbert problem formulated in the complex $k$-plane. Various Deift-Zhou contour deformations and the motivation behind them are given. By applying the nonlinear steepest descent method to the associated matrix-valued Riemann-Hilbert problem, we obtain the exact leading-order asymptotic formulas and uniform error estimates for solutions to the Cauchy problem of the Heisenberg ferromagnet equation.

We study Fredholm determinant solutions to the first two positive and the first negative members of the modified Kadomtsev-Petviashvili (mKP) hierarchy, together with the bilinear two-dimensional Toda lattice (2DTL) equation. By applying isospectral reductions, we explicitly construct Fredholm determinant solutions for the modified Korteweg-de Vries (mKdV), the modified Camassa-Holm, and the focusing nonlinear Schrödinger equations. Furthermore, we demonstrate that in the non-isospectral case, the reduction from the mKP and 2DTL equations to the non-isospectral negative KdV equation yields solutions characterized by the deformed Bessel determinant.

We study the monic recurrence coefficients of the generalized Charlier polynomials associated with the discrete weight \[ \rho_\mu(k) = \frac{\mu^k}{(k!)^2}, \qquad k\in\mathbb N_0, \qquad \mu>0. \] The coefficients $\beta_n$ and $\gamma_n$ in the recurrence relation \[ P_{n+1}(x) = (x-\beta_n)P_n(x)-\gamma_nP_{n-1}(x) \] are described through complementary moment, determinant, discrete dynamical, and continuous integrable formulations. We give a consistent normalization of the Markov function and its Jacobi continued fraction, and express the recurrence coefficients in terms of Hankel determinants, using a bordered determinant for the diagonal coefficient $\beta_n$. Starting from the original Laguerre–Freud equations, we derive the local recursive system used for computation and relate it to the Toda deformation. We also recall the known Painlevé~V representations associated with the two-parameter generalized Charlier family, in a normalization consistent with the present specialization. Moment/Hankelconstructions provide complementary moment-theoretic realizations of the recurrence coefficients, while the local Laguerre–Freud recursion provides an independent dynamical construction. Their comparison is used for numerical cross-validation. The numerical results are consistent with the leading large-degree behavior \[ \beta_n\sim n, \qquad \gamma_n\sim\mu. \]

We report on the recent results in the theory of discrete integrable equations reduced modulo a prime number $p$ based on the previous results. This article mainly concerns the $q$-discrete analogues of the Painlevé equations. We first give a quick review of the idea of `almost good reduction', which is a property that most integrable systems over a finite field $F_p$ possess when reduced from those defined over the field of $p$-adic numbers $Q_p$. This idea has been introduced in our previous papers and the discrete and q-discrete Painlevé II equations have been proved to have `almost good reduction', which can be interpreted as an arithmetic analogue of the singularity confinement test. In this article we prove that this property holds for several other $q$-discrete Painlevé equations. The theory of Diophantine integrability is reviewed and its relation to the integrability detection by the algebraic entropy is discussed.

We introduce a new method for constructing Zoll metrics that admit first integrals polynomial in the momenta. The method uses Funk's form of the metric together with an algebraic relation between the first integrals. This reduces the problem to an algebraic system from which all unknown functions can be determined explicitly. In the case of a non-trivial cubic integral, we present an explicit two-parameter family of metrics that covers all such rotationally symmetric Zoll metrics on the $2$-sphere: every such metric is isometric to a member of this family. We compare these metrics with the earlier results of Valent, Duval, and Shevchishin.

We study the Cauchy problem for the Tzitzéica equation, which is an important integrable model arising in affine differential geometry and characterizing proper affine spheres. Recently, Huang, Wang and Zhu (Math. Ann. 395, 18 (2026)) reported the long-time asymptotic result for the Tzitzéica equation for the case of the purely continuous spectrum. Inspired by their work, we conduct in-depth research on the soliton resolution conjecture and asymptotic analysis of the Tzitzéica equation. Compared with the previous work, by developing the $\bar{\partial}$-nonlinear steepest descent method, we reveal that the solution of the Tzitzéica equation exhibits three different asymptotic regions depending on $\xi:=x/t$. For the region $\xi\in(-\infty,-1]\cup[1,\infty)$, where there are no stationary phase points, we rigorously prove that the solution of the Tzitzéica equation decays algebraically to zero with error $O(t^{-1})$. In the region $\xi\in(-1,1)$, the phase function $\theta(z)$ has two stationary phase points. The corresponding asymptotic approximations can be characterized with an $N$-soliton solution as well as an interaction term between soliton solutions and the dispersion term with a residual error of order $O(t^{-3/4})$. In the transition regions $\{(x,t):1-\varepsilon<|\xi|<1\}$, the corresponding asymptotic approximation can be characterized by the $N$-soliton solution together with the interaction terms between the soliton solution and the dispersion term, with a residual error of order $O(t^{-3/4})$. Our results provide a rigorous analysis of the long-time asymptotic behavior of solutions of the Tzitzéica equation in different regions and confirm the soliton resolution conjecture for the Tzitzéica equation.

We explicitly construct the additional integrals of motion, ensuring maximal superintegrability of the discrete-time rational Ruijsenaars-Schneider model. Using them, we investigate the algebraic aspects of superintegrability in both continuous- and discrete-time settings. In particular, we determine the complete structures of the polynomial symmetry algebras associated with both the rational Ruijsenaars-Schneider model and its discretization. We demonstrate that discretization leads to a nontrivial deformation of the continuous symmetry algebra with respect to the discretization parameter, thereby extending recent analogous results from the rational Calogero-Moser system to its relativistic generalization.

We consider a hierarchy of models that can be obtained from the original Nyzhnyk system by imposing conditions on parameters, introducing potentials or pseudopotentials, performing limiting processes with respect to a scaling parameter, applying differential substitutions and interpreting parts of independent and/or dependent variables as complex or real. Therefore, among the Nyzhnyk models, one can distinguish between symmetric and asymmetric, dispersive and dispersionless, standard and modified, as well as real, complex, mixed and specific models. One can also distinguish single partial differential equations or systems of such equations, as well as linear or nonlinear Lax representations in the dispersive or dispersionless cases, respectively. For each specified model, we compute the maximal Lie invariance pseudoalgebra and, in the symmetric case, find the point- and contact-symmetry pseudogroups using the megaideal-based version of the algebraic method. It is shown that relations between these pseudoalgebras and between these pseudogroups are induced by relations between the corresponding models. Defining (finite-dimensional) subalgebras are singled out in all these pseudoalgebras in the symmetric and specific cases. Based on the established correspondences between the dispersionless and dispersive models, we completely classify one- and two-dimensional subalgebras of the maximal Lie invariance pseudoalgebra of the (dispersive symmetric potential) Nyzhnyk equation and one-dimensional subalgebras of the maximal Lie invariance pseudoalgebra of its linear Lax representation.

We obtain a large number of exact hyperbolic, trigonometric, and periodic solutions in terms of Jacobi elliptic functions as well as algebraic solutions with a power law tail of the integrable local Fokas-Lennels equation and integrable nonlocal Fokas-Lennels equation. Further, we consider a one-parameter family of generalized Fokas-Lenells equations and obtain a few of their exact solutions.

Although integrability does not universally require a continuous one-site symmetry, we rigorously prove that every jointly analytic, regular Yang-Baxter deformation of the qutrit exchange interaction necessarily retains a nontrivial, analytically varying one-site charge. Breaking this local symmetry imposes a fundamental physical constraint on approximate conservation, governed by the optimal uniform bound $\delta^3 \le C\varepsilon$ that explicitly relates the minimal one-site symmetry defect $\delta$ to the local current-conservation residual $\varepsilon$. While breaking all one-site charges strictly forbids an exact integrable completion, an optimally compensated nearest-neighbor interaction saturates this cubic limit and anomalously extends the guaranteed infinite-temperature energy-current correlation window to order $|\lambda|^{-3}$ in the perturbation strength $\lambda$. Furthermore, we reveal a fundamental resonance obstruction for intrinsic conversion perturbations that strictly prevents any exact first-order repair of a broken one-site charge on any finite ring. Nevertheless, we demonstrate that the complete eight-dimensional charge memory matrix remains thermodynamically protected and approaches the identity for timescales $t=o(|\lambda|^{-3/2})$, a robust feature of the full infinite-temperature dynamics when the thermodynamic limit is taken before weak coupling.

We study the quantized six-vertex (Q6V) model on arbitrary admissible diagrams at $q=\varepsilon$, where $\varepsilon$ is a primitive $N$th root of unity with $N$ odd. For the commuting layer transfer matrices $\mathbb{T}(y)$ with mixed boundary conditions, we establish the quantum Frobenius property \[ \mathbb{T}(y)\mathbb{T}(\varepsilon y)\cdots \mathbb{T}(\varepsilon^{N-1}y) = \mathscr{P}_N(y^N)\mathbb{I}, \] where $\mathscr{P}_N$ is a scalar Laurent polynomial obtained directly from the monomial expansion of $\mathbb{T}(y)$. Our proof exploits the underlying three-dimensional integrability, in particular the tetrahedron equations and their consequences. The quantum Frobenius relation essentially determines the joint eigenvalues of the associated Hamiltonians up to multiplicities and shows that they take a free-parafermion form, thereby extending such structures from previously known one-dimensional examples to a broad class of genuinely two-dimensional quantized six-vertex models. The construction also admits independent local gauge choices and recovers the long-range-interacting $\tau_{\; 2}$ model, which includes the original free-parafermion model as a special case. We also apply the same framework to the relativistic quantum Toda chain at odd roots of unity and obtain the free-parafermion spectra for the relevant Q6V transfer matrices with fixed boundary states.

We study the dispersionless limit of the Manakov system, the integrable two-component generalization of the nonlinear Schrödinger equation. We derive the resulting four-component genus-zero Manakov-Whitham system, characterize its hydrodynamic structure, and show that it passes the Haantjes tensor test for integrability. We show that the branch points of the spectral curve associated with plane wave solutions of the Manakov system are the local Riemann invariants of the dispersionless system. We also use the characteristic speeds to classify the baseband modulational stability/instability of the plane waves, and we study a direct linearization of the Manakov system to characterize their finite-wavenumber stability and verify agreement with the Whitham prediction in the long-wave limit. Finally, we validate the predictions by comparing them with the results of direct numerical simulations.

We investigate the equilibrium thermodynamics of a Dyson gas in connection with a set of integrable statistical mechanical observables satisfying the Toda Lattice hierarchy. We prove that in the thermodynamic limit, the integrable observables are state functions satisfying a set algebraic equations of state in closed form, obtained from direct integration of the Toda Lattice hierarchy in the continuum limit. We then explore the connection between regularity and critical behaviour of the state functions and the Dyson gas particle statistics via Monte Carlo simulations. We show that the properties of the integrable observables, such as regularity, multivaluedness, cusp singularities, carry information on the macroscopic particle statistics and its qualitative changes but with some limitations.

We compute the effective six-wave interaction in spatial and temporal Hamiltonian Dysthe equations and compare it with the complete degree-six coefficient of the deep-water Euler Hamiltonian. On every compact subset $K$ of a generic quadratic-dispersion resonance patch, exact resonances persist under a small cubic correction to the dispersion, and \[ W_6\bigl(\Phi_{\lambda,\varepsilon}(z)\bigr) =\varepsilon\frac{4BD}{aL(z)}+O_K(\varepsilon^2) \] uniformly in the resonance point $z$ and the deformation parameter $\lambda$. Here $L$ is the spectral diameter and $B,D$ are the local cubic and mean-flow coefficients. Thus $BD\ne0$ gives nonzero scattering on relatively open resonance sets. At the Craig–Guyenne–Sulem temporal coefficients, the corresponding exact Euler resonances satisfy \[ \frac{\mu}{4}W_6^{\rm E}\bigl(1+\mu\Psi_\mu(z)\bigr) =-\frac{8}{L(z)}+O_K(\mu), \] where $\mu$ is the relative bandwidth. The leading term comes from nine contractions of quartic vertices; every other connected degree-six contribution is uniformly bounded. We classify the coefficient values for which the generic six-wave coefficient vanishes identically; in the temporal case, this also gives cancellation on every non-permutation four-wave resonance. Outside these sets, the Zakharov–Schulman condition restricts every regular $C^3$ quadratic leading symbol to $\operatorname{span}\{1,k,\omega(k)\}$, the span of the symbols of wave action, momentum, and linear energy. Consequently the physical Dysthe models cannot support an inverse-scattering hierarchy with infinitely many linearly independent regular $C^3$ quadratic leading symbols. The cancellation sets contain known nonlinear Schrödinger, mixed Chen–Lee–Liu, Calogero–Moser derivative nonlinear Schrödinger, and Hirota representatives.

In this work, we develop a generalized Hamiltonian formalism for spatially nonlocal field theories whose Lagrangian densities depend explicitly on both the local field and its spatially reflected counterpart. Starting from a generalized variational principle, we derive generalized Euler-Lagrange equations and introduce a generalized functional derivative that consistently accounts for reflected-field contributions. The proposed formalism is applied to three spatially nonlocal nonlinear Schrödinger equations. For the Ablowitz-Musslimani equation, we construct, to the best of our knowledge, the first standard Lagrangian density formulated directly in terms of the complex fields and derive its complete Hamiltonian formulation. The same framework is subsequently applied to the two Lagrangian nonlocal nonlinear Schrödinger equations introduced by Velasco-Juan and Fujioka, yielding consistent Hamiltonian formulations that reproduce the corresponding generalized Euler-Lagrange equations. These results establish a unified Hamiltonian framework for a broad class of spatially nonlocal nonlinear field theories.

In this paper, we explore a variable-coefficient Fokas-Lenells system with gain/loss. Applying a nonlocal symmetry reduction, the system leads to a PT-invariant variable-coefficient Fokas-Lenells equation (vcFLE) with specific constraints on dispersion and nonlinearity coefficients. We prove the Lax integrability of the given system and construct one- and two-soliton solutions using a nonstandard Hirota bilinearization method. We show that, depending on the soliton parameters, the solutions can retain or violate the underlying symmetry. With appropriate choice of dispersion and nonlinearity, the system shows periodic-wave and soliton structures.

We study the numerical evaluation of contour integral representations for solutions of the linearized second and fourth Painlevé equations. The construction combines a nonlinear background solution, continuation of canonical Lax-pair columns, normalization matching, and quadrature over cycles with decaying branches in distinct sectors. Accuracy is assessed using a reference fundamental matrix, differential-equation residuals, and variations of monodromy data computed independently from the spectral problem. An exact identity describes the propagation of initial-basis errors. A stagewise a posteriori refinement criterion uses increments in monodromy variations to allocate numerical accuracy. The contour algorithm is compared with the Dormand–Prince method on a regular background, in a region of rapid solution growth, and through hard loss of stability followed by approximately regular oscillations. Accumulated offsets in monodromy variations and subsequent drift are measured separately. The experiments show how drift depends on spectral conditioning, quadrature accuracy, and error allocation between initial and subsequent parts of the computation. Five additional initial-data and parameter cases demonstrate the dependence of comparative performance on the nonlinear background.

In 2020 M. Mucciconi and L. Petrov introduced a long-range deformation of the quantum open non-relativistic Toda system. We prove the integrability of the deformed Toda system by constructing a $2 \times 2$ Lax operator, which produces the commutative family of differential operators containing the Hamiltonian of the deformed Toda system. Moreover, we show that the same integrable deformation exists on both classical and quantum levels and can be applied to both non-relativistic and relativistic Toda systems. For the open non-relativistic deformed Toda systems we also present an $n \times n$ Lax matrix and prove that it produces the same family of Hamiltonians. We also show how to obtain the van Diejen-type deformed Toda system. Lastly, we show that on the quantum level the algebraic Bethe ansatz technique can be applied to the deformed Toda system.

By constructing Lax operators in the loop algebra of the $\mathbb{Z}_2^2$-graded extension of the Lie superalgebra $\mathfrak{osp}(1|2)$, we derive a $\mathbb{Z}_2^2$-graded extension of the Camassa-Holm equation. The resulting equation is an integrable nonlinear PDE for a system of four $\mathbb{Z}_2^2$-graded commutative functions, each associated with a distinct $\mathbb{Z}_2^2$-degree. We further show that the $\mathbb{Z}_2^2$-Camassa-Holm equation admits a bi-Hamiltonian structure. As a consequence, it possesses infinitely many conserved quantities, including one with non-trivial $\mathbb{Z}_2^2$-degree, which are mutually in involution with respect to the $\mathbb{Z}_2^2$-graded Poisson brackets.

Many examples of integrable sigma models admit families of integrable deformations that involve coupling the undeformed “seed” theory to auxiliary fields with algebraic equations of motion. In this work, we apply this paradigm to the case where the seed theory is the lambda model, which interpolates between the Wess-Zumino-Witten model and the non-Abelian T-dual of the principal chiral model. We present an auxiliary field deformation of the lambda model and exhibit a Lax connection whose flatness is equivalent to the remaining equations of motion after imposing the algebraic field equations. We also carry out a Hamiltonian analysis of the deformed theory, verifying that it obeys the assumptions of theorems in arXiv:2605.18213 which guarantee both (i) that the model enjoys the Maillet structure that ensures involution of its conserved charges, and (ii) that the theory possesses a classical Yangian. Thus, as in other cases, the auxiliary field deformation of the lambda model changes its classical dynamics while leaving the integrable structure essentially unchanged.

We develop a geometric framework in which classical thermodynamics is reformulated as a Hamiltonian dynamical system. An explicit canonical mapping $(q,p,t,H)\leftrightarrow(V,-P,S,-T)$ identifies Maxwell's relations as the characteristic equations of the thermodynamic Poincaré–Cartan one-form, in direct parallel with Hamilton's equations of motion. Treating thermodynamic potentials as action functionals, a variational principle recovers both the Maxwell relations and the thermodynamic constraints (adiabaticity, isothermality) as conserved first integrals. For adiabatic processes in an ideal gas, this yields explicit second-order ordinary differential equations for $V(T)$ and $P(T)$, each governed by a temperature-dependent Lagrangian. Assembling the component Lagrangians into a multi-parameter Lagrangian 1-form $\mathcal{L}$, we prove that the closure condition $d\mathcal{L}=0$ holds on the solution manifold, establishing that path-independence of thermodynamic state functions is a geometric consequence of multi-time integrability rather than an independent postulate.

We revisit the Nakamura conjecture, which relates the Tomimatsu-Sato solutions of stationary axisymmetric gravity to finite Toda molecules. While the conjecture has been established partially, its general rotating sector remains an open problem. We show that the Toda determinants underlying the conjecture possess a natural weight grading. In particular, the two functions entering the Ernst potential have weights n^2 and n^2-1, and this grading extends systematically to shifted determinants labelled by partitions. In coordinates adapted to the Toda generators, each differentiation corresponds to adding one box to the associated Young diagram and increases the weight by one. The same integer n^2 also appears in the zero-order term of the Nakamura bilinear operator, revealing a compatibility between the differential equation and the determinant grading. The partition structure further explains the previously unresolved behavior of second derivatives. Repeated differentiation in one direction produces an internal sector and an external sector requiring only a one-step extension of the Wronskian hierarchy; the latter is reduced by a local three-term Pluecker relation. Thus weight grading, Young-diagram growth, Wronskian enlargement, and Pluecker reduction emerge as parts of a single determinant structure. The unit weight relation n^2 = (n^2-1) + 1 also singles out the elementary Toda seed as a natural third object, suggesting a possible route toward a genuine trilinear formulation. Although no trilinear closure is assumed here, the present construction reduces the remaining general-n Nakamura problem to definite determinant-minor identities and provides a structural framework in which such a formulation can be investigated.

We study the quantised version of the Kadomtsev-Petviashvili equation in the formulation proposed in [4]. We analyse two aspects. One is the conservation of the first higher charge $H_3$ (besides the total mass $H_0$, the momentum $H_1$ and the Hamiltonian $H_2$). We find the expression of $H_3$ and prove commutativity, provided we regularise a specific ultraviolet divergence. The result then does not depend on the choice of a regulator, under fairly general assumptions. We then study the Bethe ansatz in finite volume and obtain a general solution to the two-body problem. We have used AI supervised by us to aid the proofs.

Connection formulas for the asymptotic solutions of a system of n> 1 coupled Painleve'-II equations with symmetry-breaking parameters are written explicitly. An asymptotically exact WKB approach to these formulas relies on the quantum-mechanical independent crossing approximation for an explicitly time-dependent Schroedinger equation.

We investigate local conserved quantities in non-Hermitian bosonic lattice systems whose underlying Yang–Baxter structure remains unclear. Focusing on a one-dimensional bosonic chain with quartic interactions of creation operators, we provide a new representation of the local conserved quantities previously constructed by Sanatani and Shiraishi. Using Fourier transformation and trigonometric identities, we systematically derive these conserved quantities and show directly that they form a mutually commuting family. We further extend the construction to interactions extending beyond a single site, cubic interactions, asymmetric hopping, and an on-site potential. Our results provide a unified framework for constructing and characterizing local conserved quantities in this class of non-Hermitian many-body systems.