1. Introduction
Stochastic conflict systems constitute a cutting-edge area in modern statistical physics, characterized by emergent collective dynamics arising from competing interactions under uncertainty. The extensive modeling of these systems encompasses a wide range of physical domains, which include but are not limited to predator-prey ecosystems, chemical reaction networks, armed conflict processes, and the statistical mechanics of social impact [1–7]. Recent progress in path integral methods [8] has facilitated the first-principles derivation of switching probabilities in such conflict systems, uncovering universal scaling laws in the vicinity of critical points [9]. This physical paradigm offers a natural framework for combat dynamics to which the Wiener path integral method will be employed in this study.
Military combat systems constitute a prototypical class of driven-dissipative conflict systems where stochasticity fundamentally governs attrition dynamics. The Lanchester’s model, formalized in 1916 [10], provides a mechanical deterministic description of force decay via coupled first-order ordinary differential equations. They have been applied extensively in analyzing competitive resource-depletion processes across multiple domains such as combat simulation, market competition, enterprise management, financial gaming, virus transmission, species evolution, and network attack and defense [11–19]. While establishing quantitative foundations for operations research, these classical formulations neglect the intrinsic fluctuations which is an essential physical characteristic of real combat [20]. Existing stochastic extensions, primarily Markovian state-transition models [21–25] and stochastic differential equations [26, 27], suffer from limited applicability, poor generalizability, and the lack of a systematic mathematical foundation. They typically fail to capture the physics-informed features because of the limitations of the discrete state-space and the memoryless assumptions.
The Wiener path integral framework [28], grounded in statistical physics [29–31], may provide a mathematically rigorous framework to overcome these limitations. In this work, we address this purpose and establish a novel nonlinear Wiener path integral formalism to model stochastic combat dynamics, building upon Lanchester’s square law. This approach conceptualizes a combat system as a complex physical system of two interacting Brownian particles. It enables the computation of the probability distribution over all possible final states transitioning from a given initial state. This computed distribution is then applied to interpret the empirical 3:1 combat rule.
The so-called ‘3:1 combat rule’ persists as a controversial topic in military analysis. As per reference [32], it is formulated as an inductive rule stating that an attacker must possess an advantage of more than 3:1 on each main axis to achieve victory, which has garnered substantial attention. Nevertheless, it has been challenged on multiple fronts. There are ambiguities in the definition of proportional units, a dearth of empirical evidence, and a lack of methodological rigor in its derivation [33, 34]. Similarly, due to the absence of reliable combat data for validation, it has been regarded as a rudimentary heuristic principle [35]. In response to these issues, several analytical methodologies have been adopted, including adaptive dynamic mode modeling [33], and stochastic extensions of the Lanchester equations [22]. Notably, the stochastic model presented in reference [22] integrated probabilistic factors, thus overcoming the limitations inherent in earlier deterministic models. Given these ongoing controversies, there is an evident necessity to reevaluate the 3:1 combat rule using other mathematical frameworks.
In the next section, we will briefly review the Lanchester model. A concise introduction to the Wiener path integral is provided in section 3 to characterize interacting Brownian particles. In section 4 , Lanchester’s square law is naturally incorporated into the Wiener path integral formalism to consider the inherent battlefield fluctuations. Subsequently, in section 5 , a Wiener path integral analysis of the 3:1 combat rule is presented. We determine the temporal evolution of the attacker’s winning probability under various parameters that govern the strength of the attacker and the defender. It is found that the effectiveness of the rule, regardless of whether it benefits the attacker or the defender, is critically contingent upon the relative combat efficiencies of the opposing forces and the tolerance for attrition, which is in line with the conclusion derived from the other stochastic combat model [22]. Section 6 is for a short summary, and two appendixes are for the derivation of the diffusion matrix of Lanchester’s square law and the Rayleigh–Ritz method to solve the Wiener path integral.
2. The Lanchester’s model
The Lanchester’s model constitutes an empirical formulation describing the attrition of opposing forces. This section presents a concise introduction to the classical Lanchester’s model (for more details please refer to reference [36]). Considering an isolated system comprising red and blue armies, Lanchester’s model can be expressed as
$\begin{eqnarray}\left\{\begin{array}{rclr}\frac{{\rm{d}}{q}_{{\rm{r}}}(t)}{{\rm{d}}t} & = & {F}_{{\rm{r}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t)) & \\ \frac{{\rm{d}}{q}_{{\rm{b}}}(t)}{{\rm{d}}t} & = & {F}_{{\rm{b}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))\end{array}\right.,\end{eqnarray}$
where qr(t) and qb(t) denote the troop strengths of the red and blue forces at time t, respectively. The function Fr,b(qr(t), qb(t)) acts as the driving force governing the dynamics of both qr(t) and qb(t), where distinct functional forms correspond to different Lanchester equations. From this perspective, the Lanchester equations share a fundamental equivalence with Newton’s second law of motion in physics, i.e. $\frac{{\rm{d}}{\boldsymbol{p}}}{{\rm{d}}t}={\boldsymbol{F}}$, where the applied force equals the time derivative of momentum. This establishes the Lanchester equations as Newton’s second law for combat systems.These troop deployment patterns are denoted as the Lanchester first linear law, second linear law, and square law, applicable respectively to ancient cold-weapon combat, modern indirect-fire targeting operations, and direct-fire engagements. The ‘forces’ acting upon the red and blue sides under first linear law and second linear law are defined as 2 ) and (3 ), the solution to equation (1 ) satisfies a linear law
$\begin{eqnarray}\begin{array}{l}{F}_{{\rm{r}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))=-{\rho }_{{\rm{b}}},\\ {F}_{{\rm{b}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))=-{\rho }_{{\rm{r}}},\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{l}{F}_{{\rm{r}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))=-{\rho }_{{\rm{b}}}{q}_{{\rm{r}}}(t){q}_{{\rm{b}}}(t),\\ {F}_{{\rm{b}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))=-{\rho }_{{\rm{r}}}{q}_{{\rm{r}}}(t){q}_{{\rm{b}}}(t),\end{array}\end{eqnarray}$
where ρr and ρb are time-invariant constants independent of qr(t) and qb(t), reflecting the combat effectiveness coefficients of the red and blue sides. The ratio of the ‘force’ magnitudes acting on both sides remains constant $\begin{eqnarray}\frac{{F}_{{\rm{r}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))}{{F}_{{\rm{b}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))}=\frac{{\rho }_{{\rm{b}}}}{{\rho }_{{\rm{r}}}}.\end{eqnarray}$
Given the ‘force’ forms equations ( $\begin{eqnarray}{\rho }_{{\rm{r}}}{q}_{{\rm{r}}}(t)-{\rho }_{{\rm{b}}}{q}_{{\rm{b}}}(t)=\,\rm{constant}\,.\end{eqnarray}$
The ‘force’ exerted on both sides under the square law is defined as 1 ) yields the temporal evolution of red and blue force strengths
$\begin{eqnarray}\begin{array}{l}{F}_{{\rm{r}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))=-{\rho }_{{\rm{b}}}{q}_{{\rm{b}}}(t),\\ {F}_{{\rm{b}}}({q}_{{\rm{r}}}(t),{q}_{{\rm{b}}}(t))=-{\rho }_{{\rm{r}}}{q}_{{\rm{r}}}(t).\end{array}\end{eqnarray}$
The solution to equation ( $\begin{eqnarray}\left\{\begin{array}{rclr}{q}_{{\rm{r}}}(t) & = & {q}_{{\rm{r}}}(0)\cosh (\sqrt{{\rho }_{{\rm{r}}}{\rho }_{{\rm{b}}}}\,t)-\sqrt{\frac{{\rho }_{{\rm{b}}}}{{\rho }_{{\rm{r}}}}}{q}_{{\rm{b}}}(0)\sinh (\sqrt{{\rho }_{{\rm{r}}}{\rho }_{{\rm{b}}}}\,t), & \\ {q}_{{\rm{b}}}(t) & = & {q}_{{\rm{b}}}(0)\cosh (\sqrt{{\rho }_{{\rm{r}}}{\rho }_{{\rm{b}}}}\,t)-\sqrt{\frac{{\rho }_{{\rm{r}}}}{{\rho }_{{\rm{b}}}}}{q}_{{\rm{r}}}(0)\sinh (\sqrt{{\rho }_{{\rm{r}}}{\rho }_{{\rm{b}}}}\,t),\end{array}\right.\end{eqnarray}$
where qr(0) and qb(0) represent the initial force strengths of the red and blue sides. This solution satisfies the following square law $\begin{eqnarray}{\rho }_{{\rm{r}}}{q}_{{\rm{r}}}^{2}(t)-{\rho }_{{\rm{b}}}{q}_{{\rm{b}}}^{2}(t)=\,\rm{constant}\,.\end{eqnarray}$
The Lanchester’s model equation (1 ) describes a simplified combat scenario between two isolated systems without reinforcements. The left-hand side quantifies the rate of attrition for one force per unit time, while the right-hand side represents the attrition inflicted by the opposing force. This model, comprising the force differential equation (1 ) and its initial conditions, exemplifies mechanistic determinism. Given the combat effectiveness coefficients ρ,β and initial force strength q(0), the force strength q(t) at any time t is uniquely determined by solving the equation, as shown in solution equation (7 ). Furthermore, it embodies an extreme mechanistic deterministic proposition of ‘the winner always wins’ arising from inherent conservation laws expressed in equations (5 ) and (8 ). Here, Q1 = ρrqr − ρbqb and ${Q}_{2}={\rho }_{{\rm{r}}}{q}_{{\rm{r}}}^{2}-{\rho }_{{\rm{b}}}{q}_{{\rm{b}}}^{2}$ constitute invariant of the system. Provided the combat effectiveness coefficients remain constant during an engagement, if initial force deployments satisfy Q1,2 > 0, the red force maintains a consistent advantage throughout the conflict.
Certain analytical outcomes derived from the Lanchester equations align with established principles of warfare, thereby constituting a fundamental tool of force attrition and conflict evolution. This model represents a prevalent mathematical framework for combat simulation research and has been applied to analyze data from historical battles including Iwo Jima, Kursk, and the Ardennes [37–42], as well as the ongoing conflict in Ukraine [43].
In practice, however, such deterministic behavior stands in stark contrast with the frequently unpredictable outcomes observed in actual warfare. Factors including weather, judgment errors, personnel performance, equipment failures, intelligence acquisition, unforeseen events, and luck can profoundly influence the course and outcome of military engagements. These stochastic effects render combat an extremely complex system characterized by pervasive uncertainty, resulting in divergence between model simulations of specific historical engagements and documented outcomes, and yielding unsatisfactory results [40, 41].
To characterize and quantify the dynamic impact of such factors on combat processes, the stochastic Lanchester equation was proposed. By incorporating stochastic processes and probabilistic structures, it aims to provide a mathematically tractable framework for modeling uncertainty in warfare. For one thing, models grounded in Markov process theory introduce the Markov property into the Lanchester’s framework [21–25], enabling derivation of the system’s probability distribution through state transition probabilities. This facilitates computation of key metrics, including the force distribution at combat termination and the expected engagement duration. For another, analytical formulations based on stochastic differential equations primarily employ tools from stochastic analysis [26, 27], such as the Itô integral. These formulations introduce stochastic noise terms into the original differential equations or directly model force dynamics as Itô stochastic differential equations. Subsequently, theories of stochastic differential equations and the Itô lemma are applied to solve for the probability distribution of combat outcomes.
Although the theory of stochastic Lanchester equations has extended the original model and mitigated certain inherent limitations, it still suffers from several shortcomings, including restricted applicability, limited generalizability, and the absence of a systematic mathematical foundation. In the present study, the aforementioned disadvantages will be surmounted through the introduction of the Wiener path integral description.
3. Introduction to the Wiener path integral
The Wiener path integral method, a mathematical formalism for characterizing the stochastic trajectories of classical particles, provides distinct advantages in investigating the dynamical behavior of complex stochastic systems [28]. By comprehensively integrating over all possible random paths, this technique constructs probability distributions governing system evolution, thereby furnishing a rigorous mathematical framework for analyzing phenomena including Brownian motion and diffusion processes. Its theoretical foundation amalgamates stochastic process theory with functional integral methodology, thus enabling not only the derivation of probabilistic laws dictating system evolution but also the natural incorporation of fluctuations and non-equilibrium characteristics.
In statistical physics, the Langevin equation of multi-dimensional Brownian particles interacting with each other and moving in the presence of a linear external force reads as 9 )) and statistical regularities embedded in generalized coordinate fluctuations.
$\begin{eqnarray}\dot{{\boldsymbol{x}}}(t)+{\boldsymbol{k}}{\boldsymbol{x}}(t)={\boldsymbol{\phi }}(t).\end{eqnarray}$
The n-dimensional vector ${\boldsymbol{x}}={({x}_{1},{x}_{2},\ldots ,{x}_{n})}^{{\rm{T}}}$ is the x-coordinate of n Brownian particles, and $\dot{{\boldsymbol{x}}}$ denotes its time derivative. The n-dimensional matrix ${\boldsymbol{k}}={({k}_{ij})}_{n\times n}$ refers to the interaction coefficients. The Gaussian white noise vector ${\boldsymbol{\phi }}(t)={({\phi }_{1},{\phi }_{2},\ldots ,{\phi }_{n})}^{{\rm{T}}}$ denotes an n-dimensional random force characterized by the expected values ${\rm{E}}\left[{\phi }_{i}(t)\right]=0$ and ${\rm{E}}\left[{\phi }_{i}(t){\phi }_{j}(s)\right]={\sigma }_{ij}^{2}{\delta }_{ij}\delta \left(t-s\right)$. The diagonal matrix ${\boldsymbol{D}}=\,\rm{diag}\,({\sigma }_{11}^{2},{\sigma }_{22}^{2},\ldots ,{\sigma }_{nn}^{2})$ describes the magnitude of the intensity of diffusive motion. It is typically derived by analyzing the Langevin equation (equation (The Langevin equation is the microscopic basis for the Fokker–Planck equation and, in particular, for the diffusion equation
$\begin{eqnarray}\frac{\partial W}{\partial \tau }+\displaystyle \sum _{i}\frac{\partial ({A}_{i}W)}{\partial {x}_{i}}-\displaystyle \sum _{i}\frac{{\partial }^{2}({B}_{ii}W)}{\partial {x}_{i}^{2}}=0,\end{eqnarray}$
where W(x(τ), τ∣x(0), 0) is the transition probability density of the Brownian particle moving from the initial spatio-temporal coordinate (x(0), 0) to the terminal coordinate (x(τ), τ). The vector A = −kx represents the deterministic influence on the Brownian motion. The matrix ${\boldsymbol{B}}=\frac{1}{2}{\boldsymbol{D}}$ is associated with the effect of the noise. The equation depicts the temporal evolution of the position probability distribution of Brownian particles. Given the initial condition $\begin{eqnarray}\begin{array}{l}W({x}_{1}(t),{x}_{2}(t),t| {x}_{1}(0),{x}_{2}(0),0)\\ \quad \mathop{\longrightarrow }\limits^{t\to 0}\delta ({x}_{1}(t)-{x}_{1}(0))\delta ({x}_{2}(t)-{x}_{2}(0)),\end{array}\end{eqnarray}$
the solution of the above equation can be formulated via the Wiener path integral (for more details, please refer to chapter 1 of reference [28]) $\begin{eqnarray}\begin{array}{l}W({\boldsymbol{x}}(\tau ),\tau | {\boldsymbol{x}}(0),0)={{\rm{e}}}^{\frac{\,\rm{Tr}\,[{\boldsymbol{k}}]\tau }{2}}{\displaystyle \int }_{{ \mathcal C }\{{\boldsymbol{x}}(0),0;{\boldsymbol{x}}(\tau ),\tau \}}\\ \quad \times { \mathcal D }{\boldsymbol{x}}(t){{\rm{e}}}^{-{\displaystyle \int }_{0}^{\tau }{\rm{d}}t{\left[\dot{{\boldsymbol{x}}}+{\boldsymbol{k}}{\boldsymbol{x}}\right]}^{T}{\left[2{\boldsymbol{D}}\right]}^{-1}\left[\dot{{\boldsymbol{x}}}+{\boldsymbol{k}}{\boldsymbol{x}}\right]},\end{array}\end{eqnarray}$
where ${ \mathcal D }{\boldsymbol{x}}(t)$ is the Wiener functional integral measure, and ${ \mathcal C }\{{\boldsymbol{x}}(0),0;{\boldsymbol{x}}(\tau ),\tau \}$ is the set of all possible trajectories for the system evolution. The integrand on the exponent ${ \mathcal L }(\dot{{\boldsymbol{x}}},{\boldsymbol{x}},t)={\left[\dot{{\boldsymbol{x}}}+{\boldsymbol{k}}{\boldsymbol{x}}\right]}^{{\rm{T}}}{\left[2{\boldsymbol{D}}\right]}^{-1}\left[\dot{{\boldsymbol{x}}}+{\boldsymbol{k}}{\boldsymbol{x}}\right]$ is the Lagrangian function, and ${ \mathcal S }={\int }_{0}^{\tau }{\rm{d}}t{ \mathcal L }(\dot{{\boldsymbol{x}}},{\boldsymbol{x}},t)$ is the action functional. It postulates that the probability density for a particle transitioning from an initial position x(0) to a final position x(t) is equivalent to the summation over the contributions from all possible paths connecting these two points. The contribution of each individual path is dictated by its respective probability weight prescribed by the action functional.Recent developments in Wiener path integral research have focused on solving the integral and extending its application beyond Gaussian noise regimes. The integration of sophisticated analytical and numerical techniques include closed-form estimation, Rayleigh–Ritz approximation, quadratic action functional expansion, free-boundary estimation, and nonlinear filtering [44–49]. Extension to non-Gaussian noise environments, particularly Lévy noise [50–55] and Poisson noise [56–59], has significantly expanded the method’s applicability. These advances have enabled the method to effectively address complex system characteristics, including high dimensionality, nonlinearity, non-Gaussian noise, and fractional calculus operators, thereby substantially expanding its applicability and modeling capabilities.
The semi-classical approximation constitutes the most established and formally concise approach for solving a Wiener path integral. This method operates by fixing terminal system states and incorporating these boundary conditions into the Euler–Lagrange equations to derive analytical approximations. For low-dimensional, non-interacting systems, solutions obtained via this approach exhibit close agreement with the exact analytical solutions. However, when applied to high-dimensional, coupled complex systems, the semi-classical framework frequently encounters limitations: computational complexity scales prohibitively with system dimensionality, often precluding practical implementation or yielding non-closed-form analytical expressions.
To illustrate a two-dimensional system with given initial state (x1(0), x2(0), 0) and final state (x1(τ), x2(τ), τ), we utilize the semi-classical approximation approach to derive an analytical solution for the path that minimizes the action. The transition probability density W(x1(τ), x2(τ), τ∣x1(0), x2(0), 0) can then be computed to high precision via an expansion around this minimal-action path. Extremizing the action functional
$\begin{eqnarray*}\delta { \mathcal S }=\delta {\int }_{0}^{\tau }{\rm{d}}t{ \mathcal L }(\dot{{\boldsymbol{x}}},{\boldsymbol{x}},t),\end{eqnarray*}$
yields the Euler–Lagrange equations $\begin{eqnarray*}\left\{\begin{array}{l}\frac{\partial { \mathcal L }}{\partial {x}_{1}}-\frac{{\rm{d}}}{{\rm{d}}t}\frac{\partial { \mathcal L }}{\partial \dot{{x}_{1}}}=0\\ \frac{\partial { \mathcal L }}{\partial {x}_{2}}-\frac{{\rm{d}}}{{\rm{d}}t}\frac{\partial { \mathcal L }}{\partial \dot{{x}_{2}}}=0\end{array}\right..\end{eqnarray*}$
Solving the equations under the specified initial conditions results in the minimal-action trajectories x1(t) and x2(t). Substituting these trajectories into the transition probability density gives $\begin{eqnarray}\begin{array}{l}W({x}_{1}(\tau ),{x}_{2}(\tau ),\tau | {x}_{1}(0),{x}_{2}(0),0)\\ \,=\,N{{\rm{e}}}^{-{\displaystyle \int }_{0}^{\tau }{\rm{d}}t{ \mathcal L }(\dot{{x}_{1}}(t),\dot{{x}_{2}}(t),{x}_{1}(t),{x}_{2}(t),t)},\end{array}\end{eqnarray}$
with N the normalization constant. This solution is the leading order of the transition probability expanded around the minimal-action path.4. Wiener path integral formalism of Lanchester’s square law
To illustrate the application of the path integral approach in military operations analysis, we use the Lanchester’s square law (equation (6 )) as a fundamental example and present its corresponding path integral formulation. Such application can be also employed to other combat models. For clarity, the model is reformulated following the convention of equation (9 )
$\begin{eqnarray}\left\{\begin{array}{l}\dot{{x}_{1}}=-\beta {x}_{2}\\ \dot{{x}_{2}}=-\rho {x}_{1}\\ \end{array}\right..\end{eqnarray}$
The x1 and x2 are strength forces of two opposing sides in a combat, with the coupling coefficients ρ, β denoting their combat efficiencies accordingly.The combat system depicted by the Lanchester’s square law can be effectively compared to a pair of Brownian particles that interact with each other. This analogy regards the two combating troops as a pair of particles, with their ‘positions’ defined by their current strength levels. A troop’s attrition rate is determined by the size of its opponent, just like a force field that governs the motion of the particles. The inevitable uncertainties in warfare, such as the fog of war or command errors, play a role analogous to the random thermal noise that affects the particles. This conceptual parallel offers a solid basis for reformulating the Lanchester equations into a two-dimensional Langevin equation, a standard framework in physics describing systems subject to both deterministic and random forces. In this formulation, the advanced mathematics of the Wiener path integral may be directly applied to compute the transition probability density W(x1(τ), x2(τ), τ∣x1(0), x2(0), 0) of the strength forces. This technique enables us to analyze the entire conflict by calculating the probability of all possible paths that the battle may take, transcending a single, predetermined outcome.
Therefore, the Langevin equation and the transition probability density of the Lanchester’s square law are thus special cases of equations (9 ) and (12 ) with the matrix
$\begin{eqnarray*}{\boldsymbol{k}}=\left(\begin{array}{cc}0 & \beta \\ \rho & 0\\ \end{array}\right).\end{eqnarray*}$
The stochastic effects inherent in warfare are proportional to operational intensity, particularly manifesting as an intensifying fog of war during large-scale engagements. To achieve this mathematical characterization of combat uncertainty, the stochastic noise must be formulated as a state-dependent multiplicative term, resulting in a diffusion coefficient $\sigma \propto \sqrt{{x}_{1}{x}_{2}}$. The formal derivation is presented in the appendix Appendix A . This yields the state-dependent nonlinear diffusion matrix
$\begin{eqnarray}{\boldsymbol{D}}=\left(\begin{array}{cc}\beta {x}_{1}{x}_{2} & 0\\ 0 & \rho {x}_{1}{x}_{2}\\ \end{array}\right).\end{eqnarray}$
Then, the Lagrangian function for the Lanchester’s square law is $\begin{eqnarray}{ \mathcal L }=\frac{1}{2}\frac{{(\dot{{x}_{1}}+\beta {x}_{2})}^{2}}{\beta {x}_{1}{x}_{2}}+\frac{1}{2}\frac{{(\dot{{x}_{2}}+\rho {x}_{1})}^{2}}{\rho {x}_{1}{x}_{2}}.\end{eqnarray}$
As compared to the existing stochastic Lanchester’s models, the Wiener path integral formulation adopted in this study demonstrates some advantages for analyzing force attrition dynamics. For one thing, the path integral framework inherently accommodates high-dimensional and nonlinear interactions without reliance on state space discretization, thereby offering superior generalizability. For another thing, it provides enhanced physical intuitiveness, being grounded in the Lagrangian formalism and aligning with principles of statistical physics. This enables direct interpretation of probability evolution through dominant paths and natural derivation of fluctuation properties via generating functionals. And thirdly, by leveraging the systematic mathematical theory of path integral, the method unifies constraint handling within a continuous framework under specific combat scenarios, thereby circumventing ad hoc state-space adjustment strategies commonly required in Markovian discrete-state models.
5. Wiener path integral analysis of the 3:1 combat rule
The 3:1 rule stipulates that the attacking force should possess at least triple the combat strength of the defender to achieve victory. This principle, derived from empirical observations, represents more or less a statistical law. Consequently, the Wiener path integral formalism provides an ideal framework for its investigation, specifically enabling the evaluation of the probability distribution given by equation (13 ) using the Lagrangian function defined in equation (16 ), under the conditions prescribed by this rule.
The x1 and x2 refer to the forces of the attacker and the defender respectively. For clarity, (x10, x20) is the initial state and (x1t, x2t) is the state at some later time t. Analogously to reference [22], the attacker’s winning probability P1w is defined to be
$\begin{eqnarray}{P}_{{\rm{1w}}}={\left(1+\frac{{f}_{1}{x}_{10}}{{f}_{2}{x}_{20}}R\right)}^{-1},\end{eqnarray}$
as the critical state occurring when the defender reaches its collapse threshold (1 − f2)x20 and the attacker remains above its metastability basin. The parameters ${f}_{1}=\frac{{x}_{10}-{x}_{1t}}{{x}_{10}}$ and ${f}_{2}=\frac{{x}_{20}-{x}_{2t}}{{x}_{20}}$ are attrition rates of the attacker and the defender respectively. In the following analysis, the magnitude of their values signifies the capacity of offensive and defensive units to withstand attrition. The ratio R is $\begin{eqnarray}R=\frac{{\int }_{\left(1-{f}_{2}\right){x}_{20}}^{{x}_{20}}W\left(\left(1-{f}_{1}\right){x}_{10},{x}_{2t},t| {x}_{10},{x}_{20},0\right){\rm{d}}{x}_{2t}}{{\int }_{\left(1-{f}_{1}\right){x}_{10}}^{{x}_{10}}W\left({x}_{1t},\left(1-{f}_{2}\right){x}_{20},t| {x}_{10},{x}_{20},0\right){\rm{d}}{x}_{1t}}.\end{eqnarray}$
The ratio α = β/ρ quantifies the relative combat effectiveness between the defender and the attacker. The magnitude of α reflects the military capabilities of both belligerents—encompassing tactical proficiency, armaments, and combat resolve—while simultaneously signifying their respective roles as defender and attacker within the conflict. When opposing forces possess comparable military capabilities, the defender typically benefits from advantages in environmental familiarity, physical endurance, and concealment, resulting in superior combat effectiveness relative to the attacker. Following the approach established in reference [22], we adopt the critical threshold ratio for a draw outcome in the deterministic Lanchester’s model as the value for α under such conditions, corresponding to the well-established 3:1 rule. A draw occurs when the invariant quantity derived from equation (8 ) satisfies ${Q}_{2}=\rho {x}_{1t}^{2}-\beta {x}_{2t}^{2}=0$, which defines the critical threshold ratio
$\begin{eqnarray}{\alpha }^{* }=\frac{{x}_{10}^{2}-{x}_{1t}^{2}}{{x}_{20}^{2}-{x}_{2t}^{2}}=\frac{1-{(1-{f}_{1})}^{2}}{1-{(1-{f}_{2})}^{2}}\frac{{x}_{10}^{2}}{{x}_{20}^{2}}.\end{eqnarray}$
The x10 : x20 = 3 : 1 rule gives the values of the critical points α* to be 9.0, 9.0 and 6.12 when (f1, f2) take the values (0.9, 0.9), (0.3, 0.3) and (0.3, 0.5) respectively. Through multiple calculations with different values of x10, it was found that the winning probability P1w is independent of the numerical value of x10. The initial state (x10, x20) is simply set to (60, 20) in the calculations.Under the semi-classical approximation, the dominant contribution to the path integral originates from the minimal-action path, typically obtained by solving the Euler–Lagrange equations. However, for system with Lagrangian like equation (16 ), the Euler–Lagrange equations often form a set of coupled, nonlinear differential equations that are analytically intractable. This complexity necessitates the use of numerical approaches to estimate the classical path.
We first implement a Path Integral Monte Carlo (PIMC) method enhanced with simulated annealing. This algorithm discretizes the time evolution into 100 steps. The search for the dominant path proceeds through 1000 Monte Carlo cycles to encourage exploration of the configuration space. The system undergoes an annealing schedule with a cooling factor of 0.95 per cycle to progressively converge toward the solution. New path configurations are proposed via random perturbations with a step size of 0.1, while an infinitesimal offset of 1 × 10−5 is applied to prevent numerical divergence. While PIMC provides a reliable benchmark, its stochastic nature prevents it from yielding an explicit functional form for the path, limiting its utility for subsequent analytical work.
As a powerful alternative, we then employ the semi-analytical Rayleigh–Ritz functional optimization technique. Originally developed for Wiener path integral with constant diffusion coefficients [44], this method is extended here to Lagrangian with a nonlinear diffusion coefficient matrix. Its implementation is described in appendix Appendix B .
The Wiener path integral method yields the transition probability density W(x1t, x2t, t∣ x10, x20, 0) as its primary outcome, expressed as a function of military strength (x1t, x2t) at time t. This approach fundamentally incorporates the stochastic nature of combat dynamics, accounting for fluctuations and uncertainties that deterministic models neglect. The probability density function (PDF) depicted in Figure 1 is W(x1t, x2t, 0.2∣60, 20, 0) with combat efficiencies (ρ, β) = (1.0, 9.0). It illustrates a broad distribution of possible outcomes rather than a single trajectory, highlighting the probabilistic essence of the engagement. This result is distinct from the prediction of the Lanchester’s square law equation (7 ) which yields an absolutely definite value (x1t, x2t) = (33, 11) at time t = 0.2. The Lanchester’s model, being purely deterministic, does not capture the variability and random influences present in real combat scenarios. Interestingly, this deterministic value resides precisely within the region of maximum probability predicted by the Wiener path integral method, suggesting that while the Lanchester result is a likely outcome, it is by no means the only possible one. The Wiener method thus provides a more comprehensive and realistic framework by quantifying the range of potential states and their associated probabilities.
Figure 1. (a) The PDF calculated via Rayleigh–Ritz method. (b) The PDF calculated via Path Integral Monte Carlo method. (c) Top view of (a). (d) Top view of (b). The parameters are t = 0.2, β = 1.0, ρ = 9.0. The point marked with an asterisk (*) is the value of (x1, x2) predicted by the Lanchester’s square law ( |
The PDF in Figure 1 is computed using the Rayleigh–Ritz method and PIMC, respectively. The results exhibit a high degree of consistency. This close agreement suggests both the validity and consistency of applying the Rayleigh–Ritz technique to the present problem. Furthermore, the Rayleigh–Ritz method provides a parameterized analytical expression for the minimal-action path through a well-defined convex optimization process. Thereby it can significantly improve the computational efficiency and is employed as the principal solution technique in subsequent analyses.
The temporal dynamics of P1w, R, and the constituent ratios of R across discrete parameter sets (f1, f2, α) are plotted in the Figure 2. They are shown with time plotted on the horizontal axis and the probability or ratio value on the vertical axis. The black dashed lines are the ratio R, and the light red dotted and light blue dash-dot curves refer to the denominators and numerators of the ratio respectively. The bold red solid curves are the normalized winning probability for the attacker P1w.
Figure 2. The temporal dynamics of P1w, R, and the constituent ratios of R across discrete parameter sets (f1, f2, α). |
As illustrated in Figure 2, the proposed Wiener path integral framework reveals fundamental distinctions from deterministic Lanchester predictions, primarily through the introduction of stochastic dynamics and the newly defined attacker winning probability. Within the deterministic model, force evolution and final outcome are computed using equations (7 ) and (8 ) which prescribe fixed results for any given force ratio α. Specifically, the attacker prevails if α < α*, fails if α > α*, and a stalemate occurs exactly at the critical threshold α = α*. This binary and deterministic classification leaves no room for uncertainty or probabilistic variation; each scenario leads to one and only one outcome based solely on the initial conditions. In contrast, our stochastic analysis shifts the focus from a single predetermined trajectory to a distribution of possible evolutions, thereby enabling the evaluation of the probability of attacker victory rather than a certain yes-or-no conclusion. This approach acknowledges that real-world combat scenarios are inherently influenced by random factors—such as situational uncertainties, operational discrepancies, and external disturbances—which are not captured by classical Lanchester’s theory. By incorporating a Wiener process into the dynamics, the model captures the intrinsic randomness of engagement processes, reflecting more realistic conditions where outcomes are not preordained. Consequently, the stochastic framework demonstrates that, under the influence of randomness, all outcome states remain possible across the entire range of α, though with strongly varying likelihoods. For instance, even when α is significantly greater than α*, there remains a non-zero—albeit potentially very small—probability that the attacker may still achieve victory due to favorable stochastic fluctuations. Similarly, near the critical threshold α*, the outcomes are particularly sensitive to random effects, resulting in a smooth transition between victory and defeat probabilities rather than an abrupt switch. This continuous and probabilistic perspective not only enriches the traditional Lanchester’s model but also provides a more flexible and insightful tool for assessing combat outcomes, constituting a key innovation of our methodology.
It is evident in the figure that, when combatants possess identical tolerance for troop losses, the probability of victory for the attacking force increases with decreasing α. Conversely, when the defender maintains greater tolerance for losses, its victory probability remains largely invariant across variations in α. For parameter configurations where α = 9.0 and (f1, f2) = (0.9, 0.9) or (f1, f2) = (0.3, 0.3), engagements evolve into stalemates under these conditions. This outcome aligns with predictions derived from the deterministic Lanchester equation. Crucially, the system exhibits distinct dynamical regimes governed by the parameter α. In the regime α < α*, the winning probability decays monotonically from an elevated initial value towards approximately 50%. This behavior indicates stable favorability towards the attacker under the 3:1 force ratio. Conversely, in the regime α > α*, the probability ascends monotonically from a reduced baseline to approximately 50%, revealing the ratio’s inherent disfavorability. At criticality (α = α*), the collective fluctuation spectrum for t > 0.5 exhibits universality across parameter sets. This confirms its origin in intrinsic critical phenomena rather than parametric differences. Furthermore, under asymmetric parameters (f1, f2) = (0.3, 0.5), the system maintains defensive advantage even when α is less than the classical critical value α*, and the 3:1 force ratio remains unfavorable to the attacker. This deviation from deterministic predictions underscores the role of stochasticity and different tolerance for losses in shaping tactical outcomes.
6. Summary
In this paper, we develop a Wiener path integral formalism to incorporate stochasticity into combat dynamics based on the Lanchester’s square law. This approach enables the calculation of the probability distribution governing the evolution of force strengths, offering a continuous description of stochastic engagements that goes beyond traditional deterministic formulations. Implementing this framework to analyze the empirical 3:1 combat rule, we quantitatively evaluate the temporal evolution of the attacker’s victory probability across varying parameters. Our principal finding indicates that the applicability of the 3:1 rule is contingent upon specific parameters, primarily hinging on the relative combat effectiveness ratio α between the opposing forces and the tolerance for attrition (f1, f2).
This study establishes a theoretical framework based on stochastic dynamics and Wiener path integral method, providing a rigorous analytical toolkit for conflict system analysis. The formalism allows for the incorporation of noise-driven dynamics and uncertainty in combat outcomes, bridging a gap between classical Lanchester’s theory and real-world operational variability. Future research will extend this methodology to incorporate diverse tactical configurations of force attrition models, such as heterogeneous forces, spatially distributed engagements, and adaptive decision-making processes. Additionally, we intend to investigate significant properties, such as phase transition behaviors and correlation structures within conflict dynamics, which may offer deeper insight into the emergence of dominance and collapse in competitive systems.
Appendix A Diffusion coefficients of Lanchester’s square law
The survival status of an individual soldier in combat can be represented by a Bernoulli random variable I, where survival corresponds to I = 1, and death or loss of combat capability corresponds to I = 0. Before the battle begins, a troop has x10 soldiers, each assigned a Bernoulli variable Ii with an initial value of 1. The Bernoulli variable of a solider who is killed at time t during the battle changes to 0 and remains at 0 for all future times. Then the random strength of the troop at time t can be expressed by the Bernoulli random variable set ${\{{I}_{i}\}}_{i=1}^{{x}_{10}}$
$\begin{eqnarray}{x}_{1}(t)={I}_{1}(t)+{I}_{2}(t)+\cdots +{I}_{{x}_{10}}(t).\end{eqnarray}$
The infinitesimal increment form of equation (A1 ) is
$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{x}_{1} & = & {x}_{1}(t+{\rm{d}}t)-{x}_{1}(t)\\ & = & {\rm{d}}{I}_{1}+{\rm{d}}{I}_{2}+\cdots +{\rm{d}}{I}_{{x}_{10}},\end{array}\end{eqnarray}$
where dIj is the infinitesimal change of Ij with value 0 or −1. Taking expected value E[ · ] on both sides of the above equation yields $\begin{eqnarray}\begin{array}{rcl}{\rm{E}}[{\rm{d}}{x}_{1}] & = & {\rm{E}}[{\rm{d}}{I}_{1}]+{\rm{E}}[{\rm{d}}{I}_{2}]+\cdots +{\rm{E}}[{\rm{d}}{I}_{{x}_{10}}]\\ & = & {\rm{d}}\bar{{I}_{1}}+{\rm{d}}\bar{{I}_{2}}+\cdots +{\rm{d}}{\bar{I}}_{{x}_{10}},\end{array}\end{eqnarray}$
where ${\rm{d}}\bar{{I}_{i}}={\rm{E}}[{I}_{i}(t+{\rm{d}}t)-{I}_{i}(t)]={\rm{E}}[{I}_{i}(t+{\rm{d}}t)]-{\rm{E}}[{I}_{i}(t)]$ is the infinitesimal variation of expected value over the interval from time t to t + dt. E[dx1] can be viewed as the attrition of the troop during the interval. The Lanchester’s square law gives $\begin{eqnarray*}{\rm{E}}[{\rm{d}}{x}_{1}]=\displaystyle \sum _{i=1}^{{x}_{10}}{\rm{d}}\bar{{I}_{i}}=-\beta {x}_{2}{\rm{d}}t.\end{eqnarray*}$
Assuming that during an infinitesimal time interval dt, only the j-th soldier experiences a casualty event, while the expected values $\bar{{I}_{i}}$ of the remaining soldiers do not have time to change, it follows that for all other soldiers, ${\rm{d}}\bar{{I}_{i}}=0$, where i ≠ j. Thus, we obtain the following expression:
$\begin{eqnarray*}{\rm{d}}\bar{{I}_{j}}=-\beta {x}_{2}{\rm{d}}t.\end{eqnarray*}$
The expected value of the infinitesimal change E[dIj] is $\begin{eqnarray}{\rm{E}}[{\rm{d}}{I}_{j}]=-\beta {x}_{2}{\rm{d}}t=-P({I}_{j}(t+{\rm{d}}t)=0| {I}_{j}(t)=1),\end{eqnarray}$
where P(Ij(t + dt) = 0∣Ij(t) = 1) is the transition probability from state Ij(t) = 1 to Ij(t + dt) = 0. The variation can be derived $\begin{eqnarray*}\begin{array}{rcl}{\rm{Var}}[{\rm{d}}{I}_{j}] & = & {\rm{E}}[{\rm{d}}{I}_{j}^{2}]-{({\rm{E}}[{\rm{d}}{I}_{j}])}^{2}\\ & = & 1\cdot P({\rm{d}}{I}_{j}=1)-{({\rm{E}}[{\rm{d}}{I}_{j}])}^{2}\\ & = & P({I}_{j}(t+{\rm{d}}t)=0| {I}_{j}(t)=1)-{({\rm{E}}[{\rm{d}}{I}_{j}])}^{2}\\ & = & \beta {x}_{2}{\rm{d}}t+O({\rm{d}}{t}^{2}).\end{array}\end{eqnarray*}$
The assumption of the independent homogeneous distribution of Ij can lead to the variation of dx1 $\begin{eqnarray}{\rm{Var}}[{\rm{d}}{x}_{1}]=\beta {x}_{1}{x}_{2}{\rm{d}}t+O({\rm{d}}{t}^{2}).\end{eqnarray}$
From the Langevin equation (9 ) and the relationship established from Gaussian white noise ${\phi }_{i}(t)={\sigma }_{i}\frac{{\rm{d}}{W}_{i}}{{\rm{d}}t}$ [60], we have A5 ) with equation (A6 ) and ignoring the infinitesimal O(dt) terms, one can have 15 ).
$\begin{eqnarray*}{\rm{d}}{x}_{1}+\beta {x}_{2}{\rm{d}}t={\sigma }_{11}{\rm{d}}{W}_{1},\end{eqnarray*}$
where dW is the infinitesimal Wiener increment. Taking the variance Var[·] on both sides and noting Var[βx2dt] = O(dt2), along with the fact Var[dW1] = dt from the Itô’s lemma, we obtain $\begin{eqnarray}{\rm{Var}}[{\rm{d}}{x}_{1}]={\sigma }_{11}^{2}{\rm{d}}t+O({\rm{d}}{t}^{2}).\end{eqnarray}$
Combining equation ( $\begin{eqnarray}{\sigma }_{11}^{2}=\beta {x}_{1}{x}_{2}.\end{eqnarray}$
With a similar derivation procedure, ${\sigma }_{22}^{2}=\rho {x}_{1}{x}_{2}$ can be derived. Therefore, for the Lanchester’s square law, the diffusion coefficient matrix can be the form as equation (Appendix B The Rayleigh–Ritz functional optimization technique
The semi-analytical Rayleigh–Ritz functional optimization technique begins with determining the objective function of the optimization problem, i.e.
$\begin{eqnarray}\min { \mathcal S }={\int }_{0}^{\tau }{\rm{d}}t{ \mathcal L }(\dot{{\boldsymbol{x}}}(t),{\boldsymbol{x}}(t),t),\end{eqnarray}$
where the displacement vector is ${\boldsymbol{x}}={[{x}_{1},{x}_{2}]}^{T}$. Then, expand the path function x(t) as the approximate form $\begin{eqnarray}{\boldsymbol{x}}(t)={\boldsymbol{\psi }}(t)+{\boldsymbol{z}}{\boldsymbol{h}}(t).\end{eqnarray}$
The ψ(t) are taken the lth Hermite interpolation polynomials and its coefficient is solved to satisfy the four boundary conditions $\begin{eqnarray*}\left\{\begin{array}{l}{\psi }_{1}=\displaystyle \sum _{i=1}^{4}{a}_{1i}{t}^{i}\\ {\psi }_{2}=\displaystyle \sum _{j=1}^{4}{a}_{2j}{t}^{j}\end{array}\right.\quad .\end{eqnarray*}$
The coefficient matrix z is a 2 × l matrix, which is the optimization variables in this optimization problem. Besides, ${\boldsymbol{h}}(t)={\left[{h}_{l}(t)\right]}_{l\times 1}$ is the test function vector which satisfies the condition h(0) = h(τ) = 0. The specific form of hl(t) is taken as the followed polynomial $\begin{eqnarray*}{h}_{l}(t)={t}^{2}{(t-\tau )}^{2}{L}_{l}(t),\end{eqnarray*}$
where Ll(t) is the l-th Legendre polynomial.Subsequently, the constrained optimization is performed with the specified parameters. This process employs a gradient-based quasi-Newton method to iteratively refine the expansion coefficients z by minimizing the discretized action functional. Starting from an initial guess, the algorithm computes the numerical gradient of the action with respect to z and updates the coefficients using BFGS-based Hessian approximations to converge toward a local minimum. Regularization is applied to prevent overfitting, while adaptive numerical integration ensures accurate gradient evaluation throughout the optimization. The entire procedure, which yields both the minimal action value and the corresponding optimal path in explicit functional form, can be conveniently implemented using MATLAB’s optimization toolbox.
Analogous to the proof process in the reference [44], it follows that the Rayleigh–Ritz optimization problem for the Wiener path integral with a nonlinear diffusion coefficient matrix is a convex optimization problem. This problem necessarily possesses a global minimum, and the functional form corresponding to this optimum represents the minimal-action path functional form.
The Rayleigh–Ritz method offers several distinct advantages over the PIMC approach. Most significantly, it provides a parameterized analytical expression for the minimal-action path through a well-defined convex optimization process, whereas PIMC only generates discrete path samples without explicit functional representation. This analytical foundation makes the Rayleigh–Ritz method substantially more computationally efficient, avoiding the extensive sampling required by Monte Carlo techniques. Furthermore, the method demonstrates superior numerical stability, free from the statistical fluctuations inherent in stochastic sampling approaches. The obtained functional representation enables direct theoretical analysis using standard mathematical tools and provides enhanced extensibility to more complex scenarios, advantages that the sampling-based PIMC method fundamentally lacks.