Solution-Dependent Hirota Bilinear Systems

• Solution-dependent Hirota bilinear systems are bilinear formulations tailored to specific solution classes, reductions, or auxiliary variables rather than universally applied to all solutions.
• They enable distinct bilinear representations for phenomena like N-bright, N-dark, and breather solitons by incorporating seed-solution selection and non-autonomous features.
• These systems provide insights into integrability, Lax pair construction, and discrete reductions, offering diagnostic criteria for both exact and near-integrable behavior.

Searching arXiv for recent and foundational papers on solution-dependent Hirota bilinear systems, discrete reductions, and related integrability criteria. Solution-dependent Hirota bilinear systems are bilinear formulations in which the effective Hirota system is tied to a chosen solution, solution class, reduction, auxiliary variable set, or non-autonomous flow, rather than being a single universal bilinearization for all solutions of a nonlinear equation. In the formulation given for nonlinear evolution equations, a nonlinear equation may have many solution-dependent Hirota bilinear systems, while all solution members of a given class are associated with a single bilinear system; for the nonlinear Schrödinger equation, the $N$-bright soliton, $N$-dark soliton, and breather classes each lead to their own bilinear systems (Albazlamit et al., 21 Aug 2025). Closely related phenomena occur in hybrid discretisations of Volterra-type systems, Schwarzian reductions, extended KP hierarchies with squared eigenfunction symmetries, and tau-function recurrences whose coefficients are themselves determined by the solution (Babalic et al., 2015, Lin et al., 2013, Hay et al., 2011, Hone et al., 2017).

1. Formal setting and principal mechanisms

Hirota bilinear systems are built from Hirota’s bilinear differential operators, for example

$$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$

and often appear in the schematic form

$P(D)\{f\cdot f\}=0$

or as coupled bilinear equations for several tau functions (Bai et al., 2014, Gürses et al., 23 Nov 2025). In the fixed, classical setting, the same bilinear equation governs an entire hierarchy of soliton solutions. In the solution-dependent setting documented across the cited literature, dependence enters through seed-solution selection, through auxiliary fields that survive discretisation, through source terms such as eigenfunctions $q,r$, through reduction constraints, or through coefficients reconstructed from the tau function itself (Albazlamit et al., 21 Aug 2025, Lin et al., 2013, Hone et al., 2017).

Manifestation	Representative system	Characteristic feature
Seed-solution dependence	NLSE	Different solution classes yield distinct bilinear systems (Albazlamit et al., 21 Aug 2025)
Auxiliary-field dependence	Discrete Volterra system	Nonlinear form contains $v_n,w_n$ tied to tau data (Babalic et al., 2015)
Non-autonomous/source dependence	Extended KP hierarchy	Bilinear identities contain $q,r$ and the shift $z-t_k$ (Lin et al., 2013)
Coefficient dependence on tau variables	Birational map reductions	$Y_n$ is defined from $\tau$ and its shifts (Hone et al., 2017)
Reduction-dependent inhomogeneity	dSKdV	Bilinear coefficients cannot be simultaneously gauged away (Hay et al., 2011)

This body of work shows that “solution-dependent” does not refer to a single mechanism. It refers to a family of constructions in which the bilinear representation is adapted to the solution sector being studied.

2. Seed solutions and solution classes

A direct formulation of the concept appears in a search method for nonlinear evolution equations based on a known exact solution. The starting substitution is

$N$0

which converts the original equation into a sum of rational terms in $N$1, $N$2, and their derivatives. Integer linear combinations of these terms are then searched numerically and symbolically for exact cancellations, after which the resulting relations are rewritten in Hirota form (Albazlamit et al., 21 Aug 2025). The central conclusion is that a nonlinear evolution equation may have many solution-dependent Hirota bilinear systems, but all solution members of a given class share a single bilinear system. For the NLSE, the moving $N$3-bright soliton class yields

$N$4

whereas the dark-soliton class and the breather class produce different bilinear equations (Albazlamit et al., 21 Aug 2025).

The same dependence on the solution sector appears in integrable semi-discrete systems. For the coupled Yajima–Oikawa system, the bright and dark short-wave solitons are governed by different bilinear systems. In the bright case, $N$5, while in the dark case $N$6, the bilinear operators are modified by parameters $N$7 with $N$8, and additional algebraic constraints are required for the reduction to be compatible (Chen et al., 2015). Both sectors are solved by pfaffians obtained from reductions of Bäcklund transformations of the semi-discrete BKP hierarchy, but the bilinear equations and reduction data differ fundamentally between the two cases (Chen et al., 2015).

These results make a precise distinction between equation-level and class-level bilinearization. The equation need not possess a unique global Hirota system; the solution class may be the invariant object.

3. Discretisation, reductions, and auxiliary variables

A particularly explicit source of solution dependence arises when only part of a bilinear system is discretised. For a general two-component Volterra system, one discretisation replaces both the dispersion and auxiliary bilinear equations. A second, hybrid discretisation changes only the “dispersion” bilinear equations and leaves the auxiliary spatial bilinear equations intact. The latter yields a nonlinear discrete system in which the evolution depends explicitly on auxiliary variables $N$9 and $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$0, so the nonlinear form is solution-dependent through quantities related to the tau functions (Babalic et al., 2015). At the same time, the $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$1-soliton tau functions retain the same structure as in the continuous-time case, with the same interaction factors and modified only by the discretised dispersion relation (Babalic et al., 2015).

The same paper gives an alternative bilinearisation of the scalar Lotka–Volterra equation in which only one bilinear equation is discretised, producing a new quadrilinear or “Schwarzian-type” difference equation and an auxiliary variable $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$2 in the nonlinear form (Babalic et al., 2015). This is a concrete example of a discrete Hirota system whose nonlinear closure requires extra, solution-tied fields.

Reduction from higher-dimensional discrete systems can also force solution dependence into the coefficients. In the bilinearisation of the discrete Schwarzian KdV equation, reduction from the discrete Schwarzian KP equation yields bilinear equations

$$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$3

for $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$4, with lattice functions $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$5 and parameter $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$6 entering explicitly (Hay et al., 2011). The key point is that, unlike the parent dSKP system, these coefficients cannot be simultaneously gauged away after reduction. The bilinear equations for dSKdV are therefore inhomogeneous in a genuinely solution-dependent manner (Hay et al., 2011).

Higher-order reductions of AKNS hierarchies exhibit a related phenomenon. For AKNS($$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$7) with $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$8, the Hirota bilinear forms are inhomogeneous and require auxiliary functions; under local and nonlocal reductions, the reduced bilinear equations contain shifted or reflected arguments such as $$D_x^n D_t^m f \cdot g = (\partial_x-\partial_{x'})^n(\partial_t-\partial_{t'})^m f(x,t)g(x',t')\big|_{x'=x,\;t'=t},$$9, and the admissible solution space may become severely restricted or even trivial for some nonlocal reductions (Gürses et al., 2020). Here the dependence is not only on the equation, but on the reduction mechanism itself.

4. Hierarchies, source flows, and tau-function dependence

In hierarchy-based constructions, solution dependence often enters through additional symmetry flows or wave functions derived from the tau function. For the extended KP hierarchy with squared eigenfunction symmetry, the auxiliary flow is

$P(D)\{f\cdot f\}=0$0

and the bilinear identities involve the wave functions $P(D)\{f\cdot f\}=0$1 together with the eigenfunction pair $P(D)\{f\cdot f\}=0$2. Representative identities include

$P(D)\{f\cdot f\}=0$3

and, in the extended hierarchy,

$P(D)\{f\cdot f\}=0$4

so the bilinear structure is explicitly non-autonomous and depends on the source functions and the shifted variable $P(D)\{f\cdot f\}=0$5 (Lin et al., 2013). The resulting Hirota equations include KP equations with self-consistent sources, and the paper emphasizes that these forms are simpler than earlier formulations (Lin et al., 2013).

The extended D-Toda hierarchy provides a closely related but more structural example. Its Hirota bilinear equations are written in terms of wave functions $P(D)\{f\cdot f\}=0$6 built from the tau function, and the coefficients entering these wave functions are derivatives or shifts of $P(D)\{f\cdot f\}=0$7 (Cheng et al., 2019). Theorems stated in that work show both directions of the correspondence: every solution of the Hirota bilinear equations determines a solution of the extended D-Toda Lax hierarchy, and every Lax solution arises from a tau function satisfying those bilinear equations (Cheng et al., 2019). In this setting, solution dependence is not an ad hoc modification but the standard mode by which the hierarchy is encoded.

A further tau-function variant appears in the relation among GUE level spacing, NLS, and Painlevé IV. There the NLS bilinear system for $P(D)\{f\cdot f\}=0$8 undergoes a similarity reduction, and the reduced tau function $P(D)\{f\cdot f\}=0$9 leads to a logarithmic derivative $q,r$0 satisfying the Painlevé IV $q,r$1-form. The paper explicitly describes this as a parameter-dependent or scaling-dependent bilinear representation tied to the chosen reduction and asymptotic sector, including the Clarkson–McLeod solution (Kakei, 2015). This is not the same mechanism as source dependence, but it is another case in which the effective bilinear system is selected by the solution sector.

5. Bäcklund transformations and Lax-pair generation

Bilinear Bäcklund transformations give one of the most classical forms of solution dependence: the bilinear system relates two tau functions $q,r$2 and $q,r$3, so its very definition is a relation between solutions. For the bilinear Boussinesq equation,

$q,r$4

the bilinear Bäcklund transformation is

$q,r$5

which maps one tau function to another and serves as the starting point for deriving a Lax pair (Bai et al., 2014).

Using

$q,r$6

the bilinear Bäcklund system is converted into a nonlinear system, and a linear problem

$q,r$7

is constructed with $q,r$8. The zero-curvature condition

$q,r$9

then recovers the Boussinesq equation and confirms integrability (Bai et al., 2014). This is a paradigm case in which solution dependence is constructive rather than obstructive: the bilinear relation between two solutions is the mechanism that produces the Lax representation.

The semi-discrete Yajima–Oikawa system fits the same pattern at the hierarchy level. Its bright and dark soliton solutions are obtained from reductions of Bäcklund transformations of the semi-discrete BKP hierarchy, and the pfaffian tau functions inherit the corresponding reduction data (Chen et al., 2015). In this sense, solution-dependent bilinear systems are not peripheral to integrable theory; they are often the route by which integrable structure is extracted.

6. Integrability criteria, partial success, and recurrent misconceptions

A central misconception is that the existence of one- or two-soliton solutions is already strong evidence for integrability. For bilinear partial difference equations on a $v_n,w_n$0 stencil, the construction of one- and two-soliton solutions is possible even for non-integrable equations, whereas the existence of a generic three-soliton solution imposes severe constraints and is taken as equivalent to integrability in Hirota’s sense (Hietarinta et al., 2012). The same distinction is emphasized in the study of higher-order bilinear forms $v_n,w_n$1: many monomial forms admit three-soliton solutions only under highly restrictive parity conditions, and monomials such as $v_n,w_n$2 or $v_n,w_n$3 do not yield genuine four-soliton solutions in general (Gürses et al., 23 Nov 2025).

Another misconception is that solution dependence necessarily signals nonintegrability. This is contradicted by several integrable examples. In the two-parameter family of birational maps studied through Laurentification, the tau function satisfies bilinear recurrences that are reductions of the Hirota–Miwa equation, but the coefficient

$v_n,w_n$4

depends on the solution and initial data (Hone et al., 2017). Nevertheless, the same construction yields Poisson brackets, Lax pairs, first integrals, and Liouville integrability (Hone et al., 2017).

At the same time, solution dependence can also record a failure of full Hirota integrability. In the two-dimensional defocusing NLS flow, the single oblique soliton fits an exact Hirota bilinear formalism, but the two-soliton ansatz produces two different formulas for the interaction parameter $v_n,w_n$5, one from each bilinear equation. These coincide only in the one-dimensional limit $v_n,w_n$6, or approximately when the angle between solitons is small or the Mach number is large; outside those limits the system is only “close” to integrability (Khamis et al., 2012). This suggests that solution-dependent bilinear compatibility conditions can function as a diagnostic of near-integrable behavior as well as exact integrability.

A final recurrent point is that class-specific bilinear systems need not be arbitrary. For the NLSE, the bright, dark, and breather sectors each have stable internal bilinear descriptions, and for breathers the full coupled bilinear form is the only surviving cancellation pattern. The paper explicitly suggests that this full coupling may be a signature of integrability (Albazlamit et al., 21 Aug 2025). Within the cited literature, solution dependence therefore appears not as a defect of Hirota theory, but as a precise descriptor of how bilinear structure is organized across solution sectors, reductions, and hierarchies.