Fully Semi-Lagrangian Scheme Overview

• Fully Semi-Lagrangian scheme is a numerical method that applies characteristic-based off-grid evaluation and interpolation to the entire operator.
• It is utilized in settings such as Hamilton–Jacobi–Bellman equations, advection–diffusion systems, and coupled forward–backward problems to uniformly treat drift and diffusion.
• The scheme preserves properties like monotonicity, conservation, and stability through carefully designed dynamic programming, mass redistribution, and flux formulations.

A fully Semi-Lagrangian scheme is a numerical discretization in which the semi-Lagrangian mechanism is applied to the whole operator, or to the whole coupled system, rather than only to a transport subterm. In the literature, the phrase is problem-dependent. For Hamilton–Jacobi–Bellman equations it denotes schemes where both drift and diffusion are approximated through shifted point evaluations and monotone interpolation; for advection–diffusion systems it denotes methods that treat both advection and diffusion in semi-Lagrangian form; for mean field games it denotes fully discrete forward–backward constructions in which the backward HJB equation and the forward transport or Fokker–Planck equation are both discretized in a semi-Lagrangian spirit (Debrabant et al., 2014, Bonaventura et al., 2020, Carlini et al., 2012). The common structure is characteristic-based off-grid evaluation, interpolation or remapping, and a fully discrete update that avoids replacing a substantial part of the model by a conventional Eulerian transport or diffusion step.

1. Terminology and scope

The label “fully Semi-Lagrangian” does not have a single universal meaning. In the HJB literature, it often means that both the first-order drift and the second-order diffusion are represented through shifted evaluations rather than local finite differences. In the advection–diffusion literature, it usually means that diffusion is treated semi-Lagrangianly as well, rather than by a separate Eulerian explicit or implicit step. In coupled forward–backward systems such as MFGs, it means that both equations are discretized in a semi-Lagrangian manner and coupled at the fully discrete level (Debrabant et al., 2014, Bonaventura et al., 2017, Carlini et al., 2014).

Some papers explicitly distinguish this usage from broader characteristic-based methods that are not fully semi-Lagrangian in the narrow sense. For the BGK model, the transport part is handled by backtracking characteristics, while the collision term is a local relaxation update; the paper itself does not systematically use the phrase “fully semi-Lagrangian,” and the most accurate description is a semi-Lagrangian BGK scheme rather than a fully semi-Lagrangian reformulation of all operators (Russo et al., 2010). Likewise, the semi-Lagrangian gas-kinetic scheme for smooth flows is Eulerian finite-volume in the update of cell averages and semi-Lagrangian only in the interface flux construction (Wang et al., 2014). The FEEC scheme for incompressible flows is hybrid: the transport of a momentum 1-form is semi-Lagrangian, while viscosity, pressure, and incompressibility remain in an Eulerian finite-element weak formulation (Tonnon et al., 2023). For the Vlasov–Poisson system, the full method is an operator-split composition of fully semi-Lagrangian transport subsolvers coupled to an Eulerian LDG Poisson solve, so it is fully semi-Lagrangian only in the split-operator sense (Rossmanith et al., 2010).

Setting	“Fully” means	Representative paper
HJB equations	Drift and diffusion both via shifted evaluations and interpolation	(Debrabant et al., 2014)
Advection–diffusion or Navier–Stokes	Advection and diffusion both in SL form	(Bonaventura et al., 2017)
Forward–backward MFG systems	Backward HJB and forward transport/FP both discretized SL	(Carlini et al., 2012)

This variability in usage is itself part of the subject. A plausible implication is that “fully Semi-Lagrangian scheme” should be read as a structural description tied to the PDE class, not as a single canonical algorithm.

2. Core construction principles

The central mechanism is off-grid evolution along characteristics, followed by interpolation or remapping. In the HJB setting, the discrete operator is a dynamic programming step. A representative example is

$$S_{\rho,h}[\mu](f,i,k) := \inf_{\alpha\in\mathbb R^d} \left[ I[f](x_i-h\alpha)+\frac12 h|\alpha|^2 \right] +hF(x_i,\mu(t_k)),$$

with interpolation

$I[f](x)=\sum_i f_i\,\beta_i(x),$

where the backtracked foot $x_i-h\alpha$ is generally off-grid (Carlini et al., 2012). In the more general second-order HJB framework, drift and diffusion are both encoded in shifted evaluations,

$$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$

with the shifts chosen to satisfy moment conditions matching the continuous operator (Debrabant et al., 2014).

In advection–diffusion problems, the fully semi-Lagrangian step replaces a single departure point by a set of displaced evaluations. For the 2D vorticity equation, the update takes the form

$$\omega_i^{n+1} = \frac{1}{4} \sum_{k=1}^4 I_p[\Omega^n] \left(z_i^{n+1}+\delta_k \right), \qquad \delta_k=\pm \delta e_j,\quad \delta=\sqrt{4\nu\Delta t},$$

so the same semi-Lagrangian logic handles both advection and viscous diffusion (Bonaventura et al., 2017). In second-order advection–diffusion–reaction systems, the second-order fully semi-Lagrangian construction uses nine generalized backward characteristics,

$$z^{n+1}_{k,i} = x_i-\frac{\Delta t}{2} \left( u(x_i,t_{n+1})+u(z^{n+1}_{k,i},t_n) \right) +\sqrt{3\Delta t}\,\sigma\, e_k,$$

combined with weights $\alpha_k$, and then couples reaction through a trapezoidal term (Bonaventura et al., 2020).

A conservative variant replaces pointwise updates by fluxes over characteristic-induced departure regions. For parabolic problems in flux form, the diffusive numerical flux is

$$F_k^n = \frac12 \left( \int_{x_k}^{x_k+\delta_k} R_q[V^n](x)\,dx - \int_{x_k-\delta_k}^{x_k} R_q[V^n](x)\,dx \right), \qquad \delta_k=\sqrt{2\Delta t\,\nu(x_k,t^n)},$$

and the update is

$$v_i^{n+1} = v_i^n + \frac{1}{\Delta x} \left( F_{i+1/2}^n-F_{i-1/2}^n \right),$$

which makes the diffusive step itself semi-Lagrangian and fully conservative (Bonaventura et al., 2015).

In coupled MFG systems, the same idea is used twice: backward for the value function and forward for the density. In the first-order case, the discrete flow

$$\Phi^\varepsilon_{i,k,k+1}[\mu] := x_i-h\,Dv^\varepsilon_{\rho,h}[\mu](x_i,t_k)$$

pushes nodal masses forward, and the forward update redistributes mass through interpolation weights,

$I[f](x)=\sum_i f_i\,\beta_i(x),$0

This is the defining “fully” coupled semi-Lagrangian structure in that setting (Carlini et al., 2012).

3. Major formulations across PDE classes

For deterministic first-order MFGs, the fully discrete scheme of Carlini and Silva discretizes the Lasry–Lions first-order system with quadratic Hamiltonian $I[f](x)=\sum_i f_i\,\beta_i(x),$1. The HJB equation is treated by backward dynamic programming with off-grid backtracking and interpolation, while the continuity equation is treated by forward mass pushing along discrete characteristics induced by the gradient of the mollified discrete value function. The equilibrium is then a fixed point of the map $I[f](x)=\sum_i f_i\,\beta_i(x),$2, and Brouwer’s theorem yields existence of a discrete solution (Carlini et al., 2012). The second-order extension allows degenerate diffusion by replacing the deterministic backtracking with a weak approximation of the stochastic increments,

$I[f](x)=\sum_i f_i\,\beta_i(x),$3

and the forward Fokker–Planck step becomes a Markov-chain-type mass redistribution formula (Carlini et al., 2014).

For HJB equations more broadly, the framework of Debrabant and Jakobsen gives a unifying class of schemes for linear and fully nonlinear Bellman equations with arbitrary degenerate diffusion. The operator approximation is constructed from shifted evaluations satisfying moment conditions, and monotone interpolation is essential. In that literature, “fully semi-Lagrangian” is especially apt because both drift and diffusion are encoded through the displacements $I[f](x)=\sum_i f_i\,\beta_i(x),$4, not by separate Eulerian stencils (Debrabant et al., 2014).

Bounded-domain HJB problems add a boundary complication: SL stencils overstep the domain. Two recent lines of work address this without abandoning the fully semi-Lagrangian structure. For oblique boundary conditions, exited characteristics are replaced by obliquely projected and pushed-back points,

$I[f](x)=\sum_i f_i\,\beta_i(x),$5

together with boundary correction terms involving the oblique distance and boundary cost (Calzola et al., 2021). For Dirichlet conditions on bounded domains, the truncation lengths $I[f](x)=\sum_i f_i\,\beta_i(x),$6, the branch weights $I[f](x)=\sum_i f_i\,\beta_i(x),$7, and the direction-mixing coefficients $I[f](x)=\sum_i f_i\,\beta_i(x),$8 are adjusted so that the boundary-corrected SL stencil remains monotone and consistent (Carlini et al., 26 Mar 2025).

In fluid dynamics, the phrase commonly refers to advection and diffusion both being treated semi-Lagrangianly. In the 2D vorticity–streamfunction formulation of incompressible Navier–Stokes, advection is handled by backward characteristics and diffusion by displaced-point averaging around the departure foot. The resulting explicit update avoids diffusion linear systems; the only linear solve is the Poisson equation for the streamfunction (Bonaventura et al., 2017). A different conservative realization appears in SCOUT, where both advection and diffusion are embedded in a characteristic-based space-time balance. For nonlinear conservation laws, the advective flux contains a nonvanishing characteristic-surface term,

$I[f](x)=\sum_i f_i\,\beta_i(x),$9

and the diffusion term is incorporated by a characteristic-based Crank–Nicolson discretization that leads to a tridiagonal system (Preda et al., 3 Nov 2025).

For nonlinear diffusive conservation laws on the torus, the fully semi-Lagrangian update may combine an implicit nonlinear advection backtracking with a deterministic two-point diffusion average,

$x_i-h\alpha$0

so that both transport and diffusion are treated within the same backtracking step (Takemura, 5 Aug 2025).

4. Structural properties

Monotonicity is one of the defining analytical advantages of many fully semi-Lagrangian schemes. In the HJB setting, monotone interpolation with nonnegative basis functions implies order preservation, and convergence can then be placed in a Barles–Souganidis framework (Debrabant et al., 2014). The bounded-domain Dirichlet and oblique-boundary schemes preserve this property by reweighting truncated branches rather than introducing nonmonotone boundary corrections (Carlini et al., 26 Mar 2025, Calzola et al., 2021).

Conservation properties depend on the formulation. In MFG transport schemes, conservation follows from the partition of unity,

$x_i-h\alpha$1

so that total mass remains equal to one and positivity is preserved because the weights are nonnegative (Carlini et al., 2012). In flux-form SL diffusion, exact conservation is built into the interface-flux difference update, and the move from pointwise interpolation to cell-average fluxes is precisely what makes the method fully conservative (Bonaventura et al., 2015). SCOUT likewise preserves conservation by construction, since every update is written in telescoping interface-difference form, including the characteristic-based diffusion correction (Preda et al., 3 Nov 2025).

Positivity preservation has been especially prominent in kinetic and transport settings. The SL-DG scheme for Vlasov–Poisson preserves cell-average positivity under the stated sampling assumptions and applies a Zhang–Shu-type scaling limiter,

$x_i-h\alpha$2

so that the cell average is unchanged while nonnegative sampled values are enforced (Rossmanith et al., 2010).

Several fully semi-Lagrangian schemes derive crucial regularity from discrete semiconcavity. In first-order MFGs, after spatial mollification one obtains

$x_i-h\alpha$3

and this feeds directly into the noncrossing estimate for the transport map and, in one dimension, into the $x_i-h\alpha$4 control of the density (Carlini et al., 2012). The degenerate second-order MFG scheme uses an analogous one-dimensional one-sided Lipschitz estimate on $x_i-h\alpha$5 to obtain the same type of density control (Carlini et al., 2014).

Stability claims vary with the problem class. The fully semi-Lagrangian vorticity scheme for 2D Navier–Stokes is described as unconditionally stable and avoids linear systems beyond the Poisson solve (Bonaventura et al., 2017). The second-order fully semi-Lagrangian advection–diffusion–reaction method is designed for large time steps and avoids the large linear systems required by implicit diffusion discretizations, but its reaction coupling is Crank–Nicolson-type and therefore A-stable rather than L-stable (Bonaventura et al., 2020). In SCOUT, unconditional stability is a numerical claim supported by experiments with CFL numbers up to $x_i-h\alpha$6 (Preda et al., 3 Nov 2025).

5. Convergence, accuracy, and dimensional restrictions

Rigorous convergence results are abundant, but they are often conditional and dimension-dependent. For first-order MFGs, the reconstructed HJB solution converges locally uniformly under

$x_i-h\alpha$7

and the full coupled convergence theorem is proved only in the scalar case $x_i-h\alpha$8, under

$x_i-h\alpha$9

because the proof requires one-dimensional $$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$0 bounds on the transported density (Carlini et al., 2012). For degenerate second-order MFGs, convergence is again proved only for $$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$1, with

$$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$2

for essentially the same reason (Carlini et al., 2014).

In the HJB framework, consistency and error analysis are tied to wide stencils. For the LISL class, the truncation error is

$$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$3

and the wide-stencil scaling $$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$4 is the price of monotonicity for degenerate, non-diagonally-dominant diffusions (Debrabant et al., 2014). The same interpolation-consistency mechanism reappears in bounded domains, where the mesh condition

$$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$5

is needed in both the oblique and Dirichlet boundary constructions (Calzola et al., 2021, Carlini et al., 26 Mar 2025).

For advection–diffusion, the fully semi-Lagrangian order can be higher. The second-order method for advection–diffusion–reaction systems proves

$$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$6

so temporal second order is obtained provided interpolation error is balanced appropriately (Bonaventura et al., 2020). For one-dimensional nonlinear diffusive conservation laws with spline or Hermite interpolation of degree $$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$7, the error estimates become

$$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$8

and

$$L_k^{\alpha}[\phi](t,x):= \sum_{i=1}^M\frac{\phi(t,x+y^{\alpha,+}_{k,i}(t,x))-2\phi(t,x)+\phi(t,x+y^{\alpha,-}_{k,i}(t,x))}{2k^2},$$9

which makes explicit how interpolation order and time step must be balanced (Takemura, 5 Aug 2025).

Not all fully semi-Lagrangian schemes are high-order in time. The 2D vorticity–streamfunction method is first order overall, with local truncation error

$$\omega_i^{n+1} = \frac{1}{4} \sum_{k=1}^4 I_p[\Omega^n] \left(z_i^{n+1}+\delta_k \right), \qquad \delta_k=\pm \delta e_j,\quad \delta=\sqrt{4\nu\Delta t},$$0

and its boundary-corrected diffusion stencil may drop to $$\omega_i^{n+1} = \frac{1}{4} \sum_{k=1}^4 I_p[\Omega^n] \left(z_i^{n+1}+\delta_k \right), \qquad \delta_k=\pm \delta e_j,\quad \delta=\sqrt{4\nu\Delta t},$$1 consistency near boundaries (Bonaventura et al., 2017). Flux-form SL diffusion is likewise first order in time in the analysis presented,

$$\omega_i^{n+1} = \frac{1}{4} \sum_{k=1}^4 I_p[\Omega^n] \left(z_i^{n+1}+\delta_k \right), \qquad \delta_k=\pm \delta e_j,\quad \delta=\sqrt{4\nu\Delta t},$$2

even though higher-order spatial reconstruction is possible (Bonaventura et al., 2015).

6. Applications, limitations, and common misconceptions

Fully semi-Lagrangian schemes now appear across several research areas. In MFGs they have been used for deterministic first-order systems, degenerate second-order systems, and price-formation models where the HJ equation, transport equation, and price-clearing relation are discretized together (Carlini et al., 2012, Carlini et al., 2014, Ashrafyan et al., 2024). In PDE-constrained control and dynamic programming they provide monotone approximations for degenerate HJB equations on $$\omega_i^{n+1} = \frac{1}{4} \sum_{k=1}^4 I_p[\Omega^n] \left(z_i^{n+1}+\delta_k \right), \qquad \delta_k=\pm \delta e_j,\quad \delta=\sqrt{4\nu\Delta t},$$3 and on bounded domains with Dirichlet or oblique boundary conditions (Debrabant et al., 2014, Carlini et al., 26 Mar 2025, Calzola et al., 2021). In fluid mechanics they support large-time-step discretizations of advection–diffusion and vorticity equations without the standard explicit parabolic restriction (Bonaventura et al., 2017, Bonaventura et al., 2020). In diffusive conservation laws and kinetic transport they enable conservative or positivity-preserving updates driven by characteristic geometry rather than explicit Eulerian transport (Bonaventura et al., 2015, Rossmanith et al., 2010, Takemura, 5 Aug 2025).

A recurrent misconception is that any method with characteristic tracing is fully semi-Lagrangian. The literature does not support that simplification. The BGK scheme of Dimarco and Pareschi is semi-Lagrangian for transport and local or implicit for collision, so whether it is called fully semi-Lagrangian depends on terminology (Russo et al., 2010). The semi-Lagrangian gas-kinetic scheme for smooth flows is not fully semi-Lagrangian because the main update of the conservative variables remains Eulerian finite-volume; only the interface distribution is constructed along characteristics (Wang et al., 2014). The FEEC method for incompressible flow is semi-Lagrangian only in the transport part, while incompressibility, pressure, and viscosity stay in a variational saddle-point solve (Tonnon et al., 2023). The Vlasov–Poisson SL-DG method is fully semi-Lagrangian in each split transport stage, but not as an unsplit solver for the full coupled field-particle system (Rossmanith et al., 2010).

Another limitation is that the analytical theory often lags behind practical use. Several MFG convergence results are one-dimensional only (Carlini et al., 2012, Carlini et al., 2014). Bounded-domain HJB schemes require delicate treatment of overstepping and often only recover convergence, not general rates (Calzola et al., 2021, Carlini et al., 26 Mar 2025). Wide stencils are intrinsic to monotone second-order SL discretizations and can become costly or geometrically awkward on fine meshes (Debrabant et al., 2014). Near boundaries, consistency can degrade unless weights and displacements are rebalanced carefully (Bonaventura et al., 2017, Carlini et al., 26 Mar 2025).

The broader significance of the fully semi-Lagrangian paradigm is therefore not a single algorithmic template, but a recurring design principle: extend characteristic-based off-grid evaluation beyond pure advection, preserve monotonicity or conservation through interpolation or remapping, and retain a fully discrete formulation that remains faithful to the underlying dynamic programming, transport, or flux structure.