# FDTD Full Wave Simulations of Reconfigurable Intelligent Surfaces

E. Colella1,2, L. Bastianelli1,2, V. Mariani Primiani1,2, F. Moglie1,2
1 CNIT, Viale G.P. Usberti 181/A, 43124, Parma, Italy
2 DII, Università Politecnica delle Marche, via Brecce Bianche 12, 60131,
Ancona, Italy

###### Abstract

This paper presents the analysis of metasurfaces, here called reconfigurable
intelligent surface. The analysis is performed by numerical simulations that
implement the finite-difference time-domain method. The metasurface has been
modeled by metallic patches interconnected by varactor diodes. The
electromagnetic source consists of randomly generated plane wave. This kind of
analysis allows us to investigate the response of the metasurface when it is
hit by a random source.

††publicationid: pubid: 979-8-3503-2400-6/23/$31.00 ©2023 IEEE

##  I Introduction

Metamaterials are artificially created materials, obtained thanks to the
periodic arrangement of dielectric and metallic elements with sub-wavelength
dimensions and spacings much smaller than the wavelength [1]. These materials
display unique electromagnetic properties and have attracted growing interest
in science and technology over the past fifteen years [2, 3, 4, 5, 6]. Thanks
to metamaterials it is possible to obtain unconventional light-matter
interaction effects, such as a negative refractive index or super lenses.
Furthermore, thanks to the correct positioning of the elements constituting
the metamaterial, it is possible to manipulate the electromagnetic fields.
However, practical 3D applications of metamaterials have been hampered by
significant structural challenges [7, 8]. For this reason, interest has
focused more on 2D implementations, called metasurfaces, which are much
simpler to implement, less expensive and more manageable, being compatible
with already established planar machining techniques on silicon.
Electromagnetic metasurfaces are based on the idea of radio frequency
transmitting/receiving antenna arrays, consisting of an array of 2D resonant
elements capable of locally modifying both the phase and the amplitude of an
incident electromagnetic wave [9, 10, 11]. Furthermore, thanks to the use of
plasmon or high refractive index dielectric effects, the dimensions of the
metasurfaces can be much smaller than the wavelength. Starting from the work
developed by Yu et al. based on nano-antennas, which is capable of modifying
the wavefront of an electromagnetic wave to obtain anomalous
reflection/refraction effects [12]. This has led to the creation of various
models to increase the bandwidth and efficiency of these devices, opening the
door to a new research sector based on 2D surface optics and photonics.
Furthermore, it has been demonstrated that it is possible to use
reconfigurable dielectric or semiconductor resonators to realize smart screens
able to control the propagation of electromagnetic fields in
telecommunications. This type of screens with the possibility of automatically
reconfiguring themselves according to the environmental radio conditions are
called reconfigurable intelligent surfaces RISs [13, 14, 15]. Reconfigurable
intelligent surfaces (RIS) are surfaces made up of antennas or radiating
elements that can be dynamically modified to control and manipulate the
propagation of electromagnetic waves in a wide range of frequencies [16]. RIS
can improve the performance of wireless communication, reduce interference,
extend coverage, increase the security of wireless networks and improve the
energy efficiency of communication devices [17]. RIS represents a promising
technology for future wireless communication networks and can be used in
various scenarios, presenting a great opportunity to enhance 5G communications
in future smart radio environments. Within this paper, we evaluated the
performance of RISs in different electromagnetic conditions by computational
electrodynamics simulations to analyze their impact on the smart radio
environment.

##  II Simulation set-up

Electromagnetic simulations involve analyzing the dynamic re-configurations of
a RIS to evaluate its focusing performance. To study the electromagnetic
response of the RIS, numerical simulations were carried out using the finite
difference time domain (FDTD) method [18]. The FDTD method is a numerical
technique used to reproduce the behavior of electromagnetic waves using a full
wave approach. This technique allows to study complex geometries of different
kinds with high precision, allowing to obtain very accurate results [19, 20].
The entire FDTD code has been implemented following the standard mathematical
procedure. Our group developed a parallel home made C code which run on
supercomputer. The simulations have been performed in a working space of
$220\times 110\times 220$ for a total of $5.52\cdot 10^{6}$ cubic cells. The
cell size is 1 mm. The electromagnetic source that hits the RIS is a plane
wave randomly generated. The generation of plane wave(s) is done by a
dedicated script and then read by the FDTD code. The whole simulated domain is
divided in three domains: i) RIS domain of $120\times 10\times 120$ cells; ii)
the total field domain of $30\times 30\times 30$ cells; iii) the scattered
field domain of $20\times 20\times 20$ cells, as represented in the Figure 3.
A detailed explanation of the plane wave generation and the implementation of
separation plane between the total field and the scattered field region are
reported in [21]. To analyze the electromagnetic behavior of metasurfaces in
5G radio environment the frequency band of 0.8 GHz – 8.4 GHz with $f_{0}$=4.6
GHz has been chosen. The angle of incidence in spherical coordinates are:
$\alpha=1.57$, $\theta=0.78$ and $\phi=0.78$ for all simulations [21].

Figure 1: RIS configuration. Representation of the simulated RIS of $10\times
10$ PEC patches of 1 cm2 connected by 180 varactor diodes, vertically and
horizontally polarized. The dielectric support of $\epsilon_{r}=3$ and 1 mm of
thickness. Figure 2: Side view of the simulated RIS. On the top are reported
the PEC patches interconnected by varactor diodes. The substrate layer is in
the middle whereas on the bottom is reported the ground plane.

The whole RIS is composed by [22]: i) metal patches; ii) varactor diodes; iii)
dielectric substrate; iv) ground plane. Figure 2 reports the side view of the
simulated RIS. The substrate, i.e. the dielectric support, has been simulated
with $\epsilon_{r}=4.4$ and $\sigma=0.0025$ S/m whereas the ground plane is
simulated as perfect electric conductor (PEC).

Figure 3: FDTD domains. The innermost domain is the RIS domain, the domain in
the middle is the total field domain while the outermost domain is the
scattered domain. Figure 4: Simulation set-up. Representation of the
simulation configuration with RIS, plane wave and detecting point P.

RIS consists of $10\times 10$ resonant structures, i.e. square metal patches
modeled by PEC. Each patch is interconnected each other by 180 varactor
diodes. The varactor diodes are designed as 1 mm3 cell connecting two patches
as shown in the Figure 1, while a single patch communicates with two other
patches. The capacitance of the varactor is simulated simply by forcing the
corresponding constitutive relationship for the fields in the FDTD cell [18].
The RIS is placed in the plane $xy$. The dielectric material inside the
working domain was air and the absorbing boundary conditions for the whole
domain are the perfect matched layers (PMLs). The values of the diodes have
been 1 pF and 0.1 pF simulating a short circuit or an open circuit. The
electric and magnetic fields have been recorded in different points $P$ as
shown in the Figure 4. The temporal step was 1.5 ps for a total of 100 periods
in all simulations. In all simulations the RIS was placed on the $xz$ plane at
the center of the workspace. The simulation set-ups are shown in the Figure 4
During the analysis we also simulated the RIS as a PEC surface.

##  III Results

In the first simulation, to test the effectiveness of the separation between
the near field and the far field, we reported the electric field on a floor
consisting of $10\times 10$ points in the case without RIS at a height of 50
cells. The result is shown in the Figure 5. The reflected field values
obtained are lower than 0.2 mV/m.

Figure 5: First simulation results. Electric field on the plane of $50\times
50$ cells at high of 50 cells without RIS. Figure 6: Second simulation
results. Electric field $P(10,100,10)$ in the case of PEC surface and in the
case of RIS with diodes at 1 pF and 0.1 pF configuration. Figure 7: Third
simulation results. Electric field $P(10,100,100)$ in the case of PEC surface
and in the case of RIS with diodes at 1 pF and 0.1 pF configuration.

In that plot are reported all the probed points, namely $100$ points where the
three Cartesian components of the electric field $(E_{x},E_{y}$ and $E_{z})$
are collected. In the second case, we reported the electric field value, only
the $E_{y}$ component, collected at the point $P(10,100,10)$. In this case we
compared the electric field when the capacitors of the RIS are set to $1.0$ pF
and $0.1$ pF respectively. Moreover, in the same plot, the two capacitor
configurations are compared to a RIS simulated by a PEC surface, Figure 6.

Figure 8: Fourth simulation results. Electric field $P(210,100,210)$ in the
case of PEC surface and in the case of RIS with diodes at 1 pF and 0.1 pF
configuration. Figure 9: Field distribution of a 2D plane of the 3D domain.
All the capacitor are set to 0.1 pF.

In this point ($P(10,100,10)$) there are no difference between the diodes
configurations. We did the same comparison in terms of electromagnetic
response of the previous configurations at the point $P(10,10,100)$, as shown
in the Figure 7. We also investigate the electric field $(E_{y})$ at point
$P(210,10,210)$ with respect to the two different capacitor configurations and
by the RIS simulated as PEC surface, Figure 8. In Figs. 7 and 8 the electric
field $(E_{y})$ has some differences based on capacitor configurations while a
higher value is with the PEC surface. In Figs. 6, 7 and 8 the main time
oscillation corresponds to the central frequency, $4.6$ GHz, of the spectrum
of the exciting pulse. Figure 9 and Figure 10 report the field distributions
for a 2D plane of the 3D simulated domain. The field distribution is evaluated
over the RIS surface. In Figure 9 all capacities are set to $0.1$ pF whereas
in Figure 10 all of them are set to $1.0$ pF. In the first configuration, the
module of the electric field distribution over the RIS is higher with respect
to the case when all capacities are set to $1.0$ pF.

Figure 10: Field distribution of a 2D plane of the 3D domain. All the
capacitor are set to 1.0 pF.

##  IV Discussion and Conclusion

In this work we reported the preliminary analysis of RIS performed by FDTD
simulations. The electromagnetic source that hits the RIS consists of a random
generated plane wave. We investigated the capability of RIS to reflect and
focus the electromagnetic field in a fixed spatial point within the simulation
domain, in particular inside the scattered field domain. Currently, two RIS
configurations have been simulated, in particular when all capacitors have set
to $0.1$ or $1.0$ pF respectively. Those two configurations have also been
compared to a RIS made by a PEC surface. Simulations campaign are on progress
and next steps are: i) optimize the diodes configurations by implementing
optimization algorithm; ii) by considering a set of plane waves that act as
random source; iii) consider more capacitors configurations. FDTD simulation
is an useful and reliable tool to speed up the analysis and the design of
RISs.

## Acknowledgment

This work has been supported by EU H2020 RISE-6G project under the grant
number 101017011. We acknowledge PRACE for awarding us access to Joliot-Curie
KNL at GENCI@CEA, France.

## References

  * [1] T. J. Cui, D. Smith, and R. Liu, Metamaterials: Theory, Design, and Applications, 1st ed. New York, NY, USA: Springer, 2009.
  * [2] F Capolino. Theory and Phenomena of Metamaterials. CRC Press; 2009.
  * [3] C. Caloz, T. Itoh Electromagnetic Metamaterials, Transmission Line Theory and Microwave Applications, Wiley-IEEE Press (2005)
  * [4] K. Achouri, B. A. Khan, S. Gupta, G. Lavigne, M. A. Salem, and C. Caloz, “Synthesis of electromagnetic metasurfaces: Principles and illustrations,” EPJ Appl. Metamater., vol. 2, 2015, Art. no. 12.
  * [5] A. M. Urbas, Z. Jacob, L. Dal Negro, N. Engheta, A. Boardman, P. Egan, A. B. Khanikaev, V. Menon, M. Ferrera, N. Kinsey et al., “Roadmap on optical metamaterials” J. Opt., vol. 18, no. 9, p. 093005, 2016.
  * [6] S. Tretyakov, “Metasurfaces for general transformations of electromagnetic fields,” Phil. Trans. R. Soc. A, vol. 373, no. 2049, p. 20140362, 2015\.
  * [7] Valentine, J., Zhang S, Zentgraf T, Ulin-Avila E, Genov DA et al. Three-dimensional optical metamaterial with a negative refractive index. Nature 2008; 455, 376-379.
  * [8] Liu, N., Guo HC, Fu LW, Kaiser S, Schweizer H et al. Three-dimensional photonic metamaterials at optical frequencies. Nat. Mater. 2008; 7, 31-37.
  * [9] Chanda, D., Shigeta K, Gupta S, Cain T, Carlson A et al. Large-area flexible 3D optical negative index metamaterial formed by nanotransfer printing. Nat. Nanotechnol. 2011; 6, 402-407.
  * [10] Berry, D. C., Malech, R. G, and Kennedy, W. The reflectarray antenna. IEEE Trans. Antennas Propagat. 1963; 11, 645-651
  * [11] McGrath, D. Planar three-dimensional constrained lenses. IEEE Trans. Antennas Propagat. 1986; 34, 46-50.
  * [12] N. Yu, P. Genevet, M.A. Kats, F. Aieta, J.-P. Tetienne, F. Capasso, and Z. Gaburro, Light Propagation with Phase Discontinuities: Generalized Laws of Reflection and Refraction, Science 334 (2011), 333.
  * [13] W. Tang, M. Zheng Chen, J. Yan Dai, Y. Zeng, X. Zhao, S. Jin, Q. Cheng, and T. Jun Cui, “Wireless communications with programmable metasurface: New paradigms, opportunities, and challenges on transceiver design,” 2019, arXiv:1907.01956. [On-line]. Available: http://arxiv.org/abs/1907.01956
  * [14] E. C. Strinati, G. C. Alexandropoulos, V. Sciancalepore, M. Di Renzo, H. Wymeersch, D.-T. Phan-Huy, M. Crozzoli, R. D’Errico, E. De Car- valho, P. Popovski, P. Di Lorenzo, L. Bastianelli, M. Belouar, J. E. Mascolo, G. Gradoni, S. Phang, G. Lerosey, and B. Denis, “Wireless environment as a service enabled by reconfigurable intelligent surfaces: The RISE-6G perspective,” in 2021 Joint European Conference on Networks and Communications and 6G Summit (EuCNC/6G Summit), 2021, pp. 562–567.
  * [15] G. C. Alexandropoulos, M. Crozzoli, D.-T. Phan-Huy, K. D. Katsanos, H. Wymeersch, P. Popovski, P. Ratajczak, Y. Bénédic, M.-H. Hamon, S. H. Gonzalez, R. D’Errico, and E. C. Strinati, “Smart wireless environments enabled by riss: Deployment scenarios and two key challenges,” in 2022 Joint European Conference on Networks and Communications 6G Summit (EuCNC/6G Summit), 2022, pp. 1–6.
  * [16] Wang, D., Tan, D. Liu, L. Particle swarm optimization algorithm: an overview. Soft Comput 22,387–408 (2018). https://doi.org/10.1007/s00500-016- 2474-6
  * [17] X. Shen, Z. Zeng, X. Liu,f RIS-Assisted Network Slicing Resource Optimization Algorithm for Coex-istence of eMBB and URLLC. Electronics 2022, 11, 2575. https://doi.org/10.3390/electronics11162575
  * [18] A. Taflove, S. C. Hagness, “Computational electrodynamics: the finite-difference time-domain method“, Artech House, 2000.
  * [19] F. Moglie, “Convergence of the Reverberation chambers to the Equilibrium Analyzed with the Finite-Difference Time-Domain Algorithm”, _IEEE Trans. Electromagn. Compat._ , 46, 2004, pp. 469–476.
  * [20] F. Moglie, L. Bastianelli, V. Mariani Primiani, “Reliable Finite-Difference Time-Domain Simulations of Reverberation Chambers by Using Equivalent Volumetric Losses”, _IEEE Trans. Electromagn. Compat._ , 58, June 2016, pp. 653–660, doi:10.1109/TEMC.2016.2548520.
  * [21] F. Moglie and A. P. Pastore, “FDTD analysis of plane wave superposition to simulate susceptibility tests in reverberation chambers,” in IEEE Transactions on Electromagnetic Compatibility, vol. 48, no. 1, pp. 195-202, Feb. 2006, doi: 10.1109/TEMC.2006.870793.
  * [22] F. Costa and M. Borgese, “Electromagnetic Model of Reflective Intelligent Surfaces,” in IEEE Open Journal of the Communications Society, vol. 2, pp. 1577-1589, 2021, doi: 10.1109/OJCOMS.2021.3092217.


# Accelerating Production LLMs with Combined Token/Embedding Speculators

Davis Wertheimer, Joshua Rosenkranz, Thomas Parnell, Sahil Suneja,
Pavithra Ranganathan, Raghu Ganti, Mudhakar Srivatsa
IBM Research

###### Abstract

This technical report describes the design and training of novel speculative
decoding draft models, for accelerating the inference speeds of large language
models in a production environment. By conditioning draft predictions on both
context vectors and sampled tokens, we can train our speculators to
efficiently predict high-quality n-grams, which the base model then accepts or
rejects. This allows us to effectively predict multiple tokens per inference
forward pass, accelerating wall-clock inference speeds of highly optimized
base model implementations by a factor of 2-3x. We explore these initial
results and describe next steps for further improvements.

##  1 Introduction

As generative Large Language Models (LLMs) grow ever larger, so too do their
capabilities. However, this growth in model size in exchange for expressive
power presents a problem for deployment in production settings: as powerful as
they are, large language models incur substantial memory and computation
overhead. Well-known models such as Llama2-13B contain 13 billion parameters,
occupying roughly 24 Gb in memory using typical 16-bit weight representations.
At inference time, this substantial computational mass is devoted entirely
toward producing one single token index at a time, in the range 0 to 32k. Thus
in the extreme case, we are incurring a 24 gigabyte overhead to produce less
than 2 bytes of information!

An obvious way to rectify this imbalance would be to predict multiple tokens
at a time. Indeed, classical NLP theory has proven that even a simple 2/3-gram
language model has great predictive capability, which tells us that learned
language models should be capable of predicting more than one token at a time
with a reasonable accuracy. However, we’d also like to refrain from altering
the base models themselves, as well-known models – such as (again) Llama-13B –
are known quantities, with desirable properties that cannot be risked in a
production setting, where predictability and reliability are paramount.

One approach to squaring this circle is speculative decoding, where a smaller
draft model or speculator is trained to predict multiple tokens given a
sequence of input. These speculative tokens are produced with low cost, and
lower accuracy than the base LLM. However, we can leverage GPU parallelism
during the LLM forward pass to evaluate the output for each of these new
tokens with minimal additional overhead. Then, by comparing the outputs to the
speculated inputs, we can accept all the predicted tokens that match the
output of the base model, while rejecting all those that don’t. In this way we
can predict multiple tokens per LLM forward pass at minimal extra cost. A
deeper explanation of speculative decoding can be found in [2, 5].

This paper describes recent efforts at IBM Research to accelerate production
LLM inference using speculative decoding with novel architectures. While prior
approaches use separate LLMs as the draft model, consuming the entire sequence
of token inputs [5, 7], Medusa [2] showed that a simple multi-headed MLP
architecture can reach the same performance. This architecture uses only the
base model’s embedding vector from the most recent token, as this embedding
vector implicitly contains the contextual information derived from the
preceding input. We introduce an improved architecture based on Medusa, with
the following contributions:

  1. 1.

We show that speculator output quality can be greatly improved by conditioning
on sampled tokens, in addition to the base model context vector.

  2. 2.

We introduce an efficient two-stage training scheme, aligning our speculators
first to base model input behavior, then to output behavior.

  3. 3.

Using this speculator training pipeline, we accelerate four highly optimized
production LLMs by a factor of 2-3x.

  4. 4.

We explore the limitations of speculative decoding in a production setting,
showing that the promised speedups diminish as baseline computation and
efficiency levels increase.

  5. 5.

We outline next steps and further areas of investigation.

Code for our work is open-sourced and available at
https://github.com/foundation-model-stack/fms-fsdp, and the speculators for
our 13B-parameter base models are released on HuggingFace:
https://huggingface.co/ibm-fms.

##  2 Related Work

The original implementations of speculative decoding rely on a smaller,
independent LLM to provide drafts for the larger base model [5, 7]. While
functional, this approach carries some practical limitations. First, by
consuming the entire past token history, even small language models can
quickly balloon to large overheads. Second, the smaller models are trained on
ground-truth text, and are therefore not aligned to the behavior of the larger
base model. Thus the draft LLM could produce suffix candidates that, while
reflective of ground truth training data, are inaccurate with respect to the
base model behavior, which is the behavior we ultimately want to reproduce.
And third, if a smaller LLM trained on the same domain, with the same
tokenizer, is not available, training it from scratch might involve
unreasonable time and effort.

The recent Medusa [2] model aims to address all three of these problems by
producing drafts via a single, multi-headed MLP. This MLP takes as input the
final embedding vector for the base model’s most recent token prediction. This
allows it to access semantic context without needing the full token history,
as the base model has already put in the effort of contextualizing this token
vector, making Medusa fast at inference time for any length of prompt. Medusa
is also relatively quick and easy to train, as it is a simple model with input
that is already semantically rich.

We identify a major shortcoming in the Medusa architecture: by using only a
context vector from the base model, Medusa fails to condition on sampled token
selections, leading to poor quality n-gram predictions. EAGLE [6] and
Recurrent Drafter [9] are concurrent work, independently arriving at the same
insight. This validates our speculator design, while in contrast to those
works, we examine specifically the performance of such models in highly
optimized production settings. In particular, our base 7B model without
acceleration produces 70% more tokens per second than the equivalent 16-bit
precision baseline in the EAGLE case study (94.9 tokens/s vs 55.1 [6]), and
170% that of Recurrent Drafter (35.6 tokens/s [9]). This leads to a very
different - and more interesting - deployment and evaluation landscape.

##  3 Speculator Design and Training

In this section we describe our speculator architecture, and then discuss
efficient training schemes.

###  3.1 Speculator Design

We started with the Medusa architecture, implementing three parallel MLP
prediction heads that predict tokens $n+1,n+2,n+3$ using the embedding vector
from the base model used to predict token $n$. Looking at the predictions,
however, we quickly identified a concerning behavior: the candidate sequences
predicted by this Medusa architecture were often very poor trigrams. This is
because the Medusa MLP uses the base model embedding vector for token $n$ as
its sole input, and completely ignores the actual token sampled from the
distribution produced by that vector. This means that the predictions of
Medusa for tokens $n+1,n+2,n+3$ are not actually conditioned on token $n$ at
all, but rather on the expectation across all possible samples of token $n$,
simultaneously. Similarly, token $n+2$ is not conditioned on token $n+1$, but
instead on the expectation across all values for $n+1$, and similarly for
token $n+3$. This leads to counterintuitive behavior.

As a simple contrived example, imagine that the base model has produced an
input ending with, “follow these steps:$\backslash$n”. Let the three most
likely continuations of this input be (without too much suspension of
disbelief required):

  * •

“Step 1:”

  * •

“$\backslash$n Step 1:”

  * •

“$\backslash$n $\backslash$n Step 1:”

Moreover, in each of these continuations, ”Step” is highly likely (as we know
it must appear), while ”$\backslash$n” is likely to have lower confidence. In
this scenario, the Medusa prediction for the next three tokens is then “Step
Step Step”, as those are clearly the three most likely tokens at each position
$n+1,n+2,n+3$. However, this is also clearly a poor predictor of base model
behavior.

To alleviate this problem, we modify the Medusa architecture into a multi-
stage MLP, where each stage predicts a single token in the draft sequence.
Each stage of the MLP takes as input both a state vector and one (or more)
embeddings for the sampled token(s) from the prior stage (the base model can
be considered stage 0). Each input pair (state plus sampled token embedding)
is projected and passed through a LayerNorm/GeLU activation, forming a new
state vector. This state vector is used to predict the next set of draft
tokens, which, with the new state vector, acts as the inputs for the next
stage of prediction. The model thus emits an entire tree of candidates, and we
allow the user to select the top-k most likely predictions from that tree to
use in practice. The 3-stage model architecture is visualized in Fig. 1.

Figure 1: A simple architecture diagram for a 3-headed multi-stage MLP
speculator. Z is the latest state vector from the base model, while T[n] is
the sampled token at time $t+n$.

At inference time, we implement multi-candidate decoding [8] and Tree
Attention [2, 7] to evaluate the chosen $k$ candidates simultaneously in
parallel. Rather than compare input tokens to output tokens for a single
candidate sequence, and taking only the correct tokens, we instead perform
this scoring process for all candidates simultaneously, and keep the tokens
from only the best-scoring candidate. This is a simple way to increase
speculator accuracy (and thus the number of generated tokens per step), at the
cost of increased parallelism. When we report tokens per step results,
therefore, we will also report the number of candidates being evaluated in
parallel.

By conditioning on prior tokens as well as the base model’s internal state,
our speculator is able to produce more realistic n-gram drafts. At the same
time, it’s also more accurate at predicting base model behavior. For
Llama2-13B, we were able to accelerate the above three-headed sandbox Medusa
architecture from 2.28 tokens per step, to 2.63 tokens per step (best of 10
candidates), simply by conditioning each MLP head on the previous sampled
tokens, in addition to the base model state.

Note that Fig. 1 could just as easily depict an RNN, and indeed we could
implement our multi-stage MLP as an RNN by simply tying the weights across the
different stages. While this weight sharing would shrink the parameter count
of our speculator, it does not make the speculator any faster at inference
time. Therefore we implement our speculator without any weight sharing for
maximum adaptability. The resulting parameter counts thus appear large
(800M-1.1B for a 7B-13B scale base LLM). But despite having $\sim 1/10$th the
parameters of the base LLM, our speculators execute in under $1/30$th of the
time, due to their shallow depth and simple architecture.

###  3.2 Speculator Training

To train our speculators, we leverage the efficient Foundation Model Stack
distributed training platform introduced in [1]. We adapt the Llama training
script to instead freeze the base Llama model and instead train the speculator
on top. We keep the same dataloaders, checkpoint utilities, and other
convenient pipeline components.

Since the speculator is trained from scratch, we have the opportunity to fully
align it to the base model: instead of predicting ground truth text, the
speculator should train directly against the generated base model output.
Unfortunately, text generation with Llama is not efficiently parallelizable –
the whole reason we’re training a speculator in the first place. We therefore
adopt a 2-stage training paradigm. In the first stage, we train on small
batches with long context lengths (4k tokens). The base model generates
embedding vectors for each token in parallel, and the speculator uses these
vectors, plus the ground truth tokens, to predict subsequent ground truth
tokens. This is a similar task to standard causal language model training,
only we’re predicting multiple tokens at each position (one for each stage of
the speculator MLP).

In stage 2, we switch to large batches with short sequence lengths (256
tokens), generated from the base model itself. By only generating a short
sequence of tokens, we can continue to train reasonably efficiently in
parallel, while fine-tuning the speculator to explicitly match base model
outputs. We did not observe any consistent overfitting to short sequence
lengths – in fact, in some cases, the speculator became more accurate on
longer inputs. In our setup, a roughly 5:2 ratio of steps for stage 1 vs stage
2 translates to a roughly 1:1 ratio in wall clock time, so we use this ratio
for convenience. The speculator efficiently trains on proxy ground-truth text
data for half the time, and then spends the other half on less efficient model
alignment.

By the end of stage 2 training, we are able to accelerate a three-headed
Llama2-13B speculator from 2.41 tokens per step at the end of stage 1, to 2.63
(best of 10 candidates). The training curves for this run across both training
stages are provided in Fig. 2.

Figure 2: Per-head training loss curves for Llama2-13B speculator training,
stages 1 (left) and 2 (right). Loss values jump downward at the start of stage
2 as the task becomes inherently easier: rather than using base model behavior
to predict ground truth text, we are using base model behavior to predict
other base model behavior (future tokens).

##  4 Results

We train four rounds of speculator models, and discuss the results for each.
The first is a proof-of-concept 3-headed speculator for Llama2-7B, which we
use to illustrate both the successes and limitations of our approach in
production settings. The second is a similar 3-headed speculator for
Llama2-13B, which we deploy to an internal LLM server. The third is a 7-headed
speculator for Codellama-13B-instruct, which achieves much better speedup due
to the more predictable nature of code. Finally, we train a 5-headed
speculator for an internal Granite-20B code model, demonstrating further
scalability with model size.

All speculators are trained on 32 NVIDIA A100 gpus, with latent dimensionality
4096, 15k steps for stage 1 training, and 6k for stage 2. For natural language
models, stage 1 batches contain 1M tokens, and stage 2 batches consist of 786M
tokens generated from 197M tokens’ worth of truncated training documents.
Because the code models are more computationally involved (7 or 5 heads vs 3)
we shrink batch sizes to 786M tokens for stage 1, and 524M/131M
generated/ground truth tokens for stage 2. All layers are initialized to
standard deviation $4096^{-1/2}$, and state/embedding projections are weighted
such that the base model’s state vector input accounts for $50\%$ of the
variance of the state at the final head, in expectation. Learning rate is set
to $1e-3$ for stage 1, with warmup over the first $5\%$ of training and cosine
decay to $1e-4$. Stage 2 follows the same LR scheme, but with a max LR of
$1e-4$ and annealing to $1e-5$, as we found that higher values caused loss to
spike badly during the transition from stage 1 to 2.

###  4.1 Llama2-7B-chat

In Table 1 we track the per-token latency of greedy inference with a Llama2-7B
model and our trained speculator on a single A100, in fp16 precision, as a
function of batch size ($b$), prompt length ($p$), and number of simultaneous
candidates ($k$). The speculator generates a tree with 6 branches at the first
head, 3 at the second, and 2 at the third, and we prune this tree of
$6*3*2=36$ candidates down to the desired $k$ base on confidence. For
reference, we include non-speculative inference under the heading $k=0$, and
the average tokens generated per step ($\tau$) across all batch sizes.
Latencies are calculated over 256 trials of exactly 100 generated tokens.
Prompts are Webhose documents truncated to the desired prompt length.

ms/tok (64) | 0 | 1 | 2 | 4 | 8 | 16 | 32
---|---|---|---|---|---|---|---
1 | 10.54 | 6.23 | 5.76 | 5.50 | 5.55 | 6.00 | 7.21
2 | 10.66 | 6.32 | 6.01 | 5.91 | 6.32 | 7.62 | 11.15
4 | 10.78 | 6.78 | 6.67 | 6.95 | 8.73 | 11.90 | 18.31
$\tau$ | 1.00 | 2.21 | 2.44 | 2.67 | 2.84 | 2.96 | 3.03

ms/tok (2048) | 0 | 1 | 2 | 4 | 8 | 16 | 32
---|---|---|---|---|---|---|---
1 | 12.87 | 9.11 | 8.89 | 8.91 | 9.68 | 11.51 | 14.78
2 | 14.98 | 10.78 | 10.94 | 11.56 | 13.24 | 16.83 | 24.96
4 | 19.62 | 14.24 | 15.63 | 16.78 | 20.97 | 29.05 | 44.88
$\tau$ | 1.00 | 2.01 | 2.21 | 2.42 | 2.61 | 2.74 | 2.83

Table 1: Iterative latency (milliseconds per token) for Llama2-7B. Rows
indicate batch size $b$, columns indicate number of parallel candidates $k$.
$k=0$ indicates non-speculative baseline, and logical speedup (tokens per
step) is given as $\tau$. Prompt length $p$ is 64 (left) vs 2048 (right). As
baseline computational load increases, speculative decoding provides less
improvement.

From the table we can make a few notable observations. First, our logical
speedups ($\tau$, increasing asymptotically with $k$) do translate to wall-
clock speedups, but not one to one. This makes sense in a production setting,
as our inference implementation is already highly optimized and torch compiled
[4]. Our baseline $b=1,k=0$ latency of 10.54 milliseconds per token translates
to 94.9 tokens per second for Llama2-7B, compared to EAGLE’s gpt-fast
implementation baseline of 55.1 tokens per second [6]. Note that speculative
decoding relies on unutilized GPU bandwidth to parallelize the evaluation of
multiple candidates and tokens during the base model forward pass. If GPU
throughput is already maximized, that leaves less resources to dedicate toward
speculative decoding based acceleration. Here we observe just under 2x speedup
for prompt length 64, batch size 1, and 4 candidates per step. This is less
than our measured tokens per step (2.67), but still a massive improvement in a
production setting!

Second, and consistent with the above intuition on limited GPU bandwidth: as
prompt length and batch size increase, the speedup from speculative decoding
erodes. For large prompts and batch sizes, speculative decoding can be over
twice as slow as non-speculative decoding, despite generating almost three
tokens per single forward pass. This is less surprising when we consider that
for $b=4,k=32$ and 4 tokens per candidate (1 ground truth token from the base
model plus 3 speculated suffix tokens), we’re running 512 tokens per step. In
terms of raw FLOPs, the equivalent non-speculative workload is a batch size of
512, far beyond the capabilities of an A100 gpu, which we observed failing at
$b\approx 128$ due to memory limits. Thus while a promising technique for
inference acceleration, the wall-clock speedup from speculative decoding
scales inversely to both the computational workload and the degree of pre-
existing optimization of the base model.

Finally, because of this falloff in performance, the optimal number of
parallel candidates for a given batch size is usually quite small in our
setup: less than 8 in all cases for prompt length 64, and less than 4 when
prompt length is 2048. This is in contrast to prior work in speculative
decoding, which usually features larger sets of candidates. For example,
Lookahead Decoding [3] obtains their best results for $k=15$, while [8]
publishes results for $k=32$. In our case, inflating the input tensor to this
size hurts more than helps, and also calls into question the contribution of
tree attention in our implementation, as there are likely to be very few
redundant tokens when $k$ is only 4 or less.

In sum, the Llama2-7B speculator results indicate that speculative decoding
produces significant material speedups in production settings. However, the
practitioner must carefully account for computational load. The choice of
hyperparameters and even algorithms may require dynamic adaptation to
workloads at runtime.

###  4.2 Llama2-13B-chat

Moving beyond proof of concept in Llama2-7B, we train a new speculator, using
the same hyperparameters, for Llama2-13B. Speculator architecture remains the
same, except that the initial state vector input now has dimensionality $5120$
instead of $4096$, as the base model width has increased. We then deployed the
accelerated base model to an internal IBM model inference server. Benchmarks
of latency and throughput for $k=[1,5]$ and varying $b$ are shown in Fig. 3.
We observe a similar 2x reduction in latency for $b=1$ and $k=5$. However,
consistent with the observed behavior for 7B, this improvement gradually
disappears for larger $b$, until speculative decoding with $k=5$ becomes
significantly slower than the baseline.

Figure 3: Throughput (x-axis) vs iterative token latency (y-axis) for
Llama2-13B with speculative decoding on an inference server. Numbers indicate
concurrent users. Top: artificial homogeneous workloads (49 tokens in, 100
tokens out). Bottom: heterogeneous workloads matched to historical logs.
Speculators scale better in the heterogeneous case (but at lower throughput)
because batch size is upper-bounded by concurrent users $b$, but in practice
is lower in expectation.

Fortunately, we can mitigate this behavior by fixing $k=1$, trading off a
slightly higher latency for $b<4$ for much better behavior at $b>4$. In the
extreme case, we can also disable speculative decoding dynamically for $b$
above a given threshold, such as 32. In this way, we can attain $>2x$ speedups
in the best case, and perfectly reproduce baseline behavior in the worst case.
These results suggest that in large-scale production contexts, a well-trained
speculator may be most performant with $k=1$, and no special multi-candidate
handling, as the marginal improvement in recall is not worth incurring the
overhead of evaluating extra candidates.

###  4.3 Codellama-13B-instruct

Shifting to a different data domain, we train a speculator for
Codellama-13B-instruct. Because code is more formalized, repetitive, and
therefore predictable than natural language, we increase the number of
stages/heads in our speculator from 3 to 7 (raising input length from 4 to 8)
to capitalize on this predictability. Performance benchmarks are given in
Table 2 for prompt lengths 64 and 512.

ms/tok (64) | 0 | 1 | 2 | 4 | 8 | 16
---|---|---|---|---|---|---
1 | 18.88 | 6.3 | 5.99 | 5.67 | 5.76 | 6.96
2 | 19.34 | 9.01 | 5.51 | 5.64 | 7.22 | 8.71
4 | 19.60 | 5.55 | 5.94 | 7.86 | 10.08 | 13.99
$\tau$ | 1.00 | 5.02 | 5.31 | 5.56 | 5.76 | 5.86

ms/tok (512) | 0 | 1 | 2 | 4 | 8 | 16
---|---|---|---|---|---|---
1 | 20.12 | 6.46 | 6.44 | 10.55 | 6.60 | 7.98
2 | 21.14 | 6.25 | 6.19 | 6.90 | 9.03 | 10.93
4 | 22.82 | 6.71 | 7.55 | 10.16 | 13.43 | 17.91
$\tau$ | 1.00 | 5.49 | 5.73 | 5.92 | 6.10 | 6.21

Table 2: Iterative latency (milliseconds per token) for
Codellama-13B-instruct. Rows indicate batch size, columns indicate number of
parallel candidates. Candidates 0 indicates non-speculative baseline, and
logical speedup (tokens per step) is given as $\tau$. Prompt length is 64
(left) vs 512 (right). As baseline computational load increases, speculative
decoding provides less improvement.

Here we observe that the increased heads are suitable for code workloads. By
doubling the maximum acceptable n-gram length, we roughly double the mean
accepted n-gram length $\tau$. This $>5x$ logical speedup produces a $>3x$
wall-clock speedup, allowing us to operate a 13B parameter model at fp16
precision with a throughput of 181.5 tokens per second ($5.51$ ms/token at
$b=2,k=2,p=64$).

Beyond the notably improved performance, results in Table 2 largely mirror
those in Table 1. Improvement falls off as batch size and prompt length
increase, shrinking the optimal $k$, which remains small in general. A
discrepancy is the fact that speculator accuracy, as indicated by $\tau$,
actually increases slightly with prompt length, as opposed to decreasing
slightly in Table 1.

###  4.4 Granite-20B

Finally, we train a speculator for a larger, internal code model, Granite-20B.
Due to memory constraints at training time, we reduce the number of speculator
heads in our code speculator from 7 to 5, and otherwise train the same as for
Codellama-13B. Similar to our Llama2-13B speculator, we then deploy this model
to our internal server and benchmark latency and throughput across batch
sizes, with results given in Fig. 4.

Figure 4: Throughput (x-axis) vs iterative token latency (y-axis) for
Granite-20B with speculative decoding on an inference server. Numbers indicate
concurrent users. Top: artificial homogeneous workloads (50 tokens in, 100
tokens out). Bottom: heterogeneous workloads matched to historical logs.
Speculators scale better in the heterogeneous case (but at lower throughput)
because batch size is upper-bounded by concurrent users $b$, but in practice
is lower in expectation.

Behavior mirrors prior results, despite the larger model size. Consistent with
Codellama-13B, we see a roughly $3x$ speedup over the non-speculative baseline
in the best case. Similar to Llama2-13B, we find that fixing $k=1$ buys us
much better scaling behavior, at the cost of a slight latency drop for very
small batch sizes. Heterogeneous workloads scale better than homogeneous ones,
though with lower throughput, due to smaller batch sizes in practice
(simulated users no longer artificially synchronize their calls). And by
disabling speculative decoding for large batch sizes, we can achieve the best
of both worlds.

##  5 Discussion

While achieving a 2-3x speedup in a production setting is significant, the
approach discussed so far represents a version 1 first attempt. A number of
avenues are available to explore for further improving the efficiency of our
speculator models.

First, our natural language speculators are constrained by their limited
number of heads. Because we only speculate 3 tokens at a time, we are capped
at 4 tokens per step (1 ground truth prediction from the base model, plus up
to 3 accepted speculated tokens). Increasing the number of heads to 5, as in
[6] and [9], is a simple and easy way to improve the performance of our
models. The same argument may also be applied to our 5-headed Granite-20B code
model, though the fact we attain the same wall-clock latency improvement as
for our 7-headed Codellama-13B speculator, suggests that heads 6 and 7 may be
redundant. Further study is needed.

Second, we declined to tie any weights in our speculators, either between
heads/stages, or between the speculator and the base model. This is because
doing so provides no wall-clock speedup. However, decreasing the parameter
count via weight tying can be useful in terms of shrinking GPU memory
overhead, especially as we increase the number of heads/stages in our
speculators. It might also lead to improved convergence, for example by
initializing the embedding layers or shared layer to the semantically rich
values of the base model’s embedding layer. Working in the same latent space
as the base model also opens opportunities for introducing useful auxiliary
losses, as in [6].

At inference time, the optimal $k$ depends on $b$, and so we’d like to swap
dynamically to the optimal value for each incoming inference request, step to
step, rather than simply toggling on/off between default values of
$k=(0,1,5)$. However, the optimal scaling plan will be model- and hardware-
specific, and thus, though desirable for an efficient speculative decoding
implementation, unlikely to yield many broader insights.

##  6 Conclusion

In this technical report we introduced a novel improvement on the Medusa
speculative decoding architecture, improving predicted token sequences by
conditioning on prior sampled tokens, in addition to the base model’s current
state vector. Using a two-stage training process, we produce trained
accelerators for Llama2-7B, Llama2-13B, Codellama-13B, and an internal
Granite-20B model, and discuss performance impacts on a highly optimized
production-level baseline. In particular, speculative decoding delivers its
promised speedups in optimal scenarios, and we are able to run blazingly fast
at fp16 precision. However, the improvement degrades for larger workloads, as
the available window for handling increased parallelism on the GPU shrinks.
Thus, delivering reliable acceleration via speculative decoding in production
settings requires thoughtful – and careful – design and effort.

##  7 Acknowledgments

We thank Supriyo Chakraborty, for insightful design discussions in the early
stages of the project, and Antoni Viros i Martin, Linsong Chu, and Brian
Vaughan, for valuable consultations on code.

## References

  * [1] Team PyTorch at IBM and Team PyTorch at Meta.  Maximizing training throughput using pytorch fsdp.  PyTorch Blog, 03 2024.
  * [2] Tianle Cai, Yuhong Li, Zhengyang Geng, Hongwu Peng, Jason D. Lee, Deming Chen, and Tri Dao.  Medusa: Simple llm inference acceleration framework with multiple decoding heads.  arXiv preprint arXiv: 2401.10774, 2024.
  * [3] Yichao Fu, Peter Bailis, Ion Stoica, and Hao Zhang.  Break the sequential dependency of llm inference using lookahead decoding, 2024.
  * [4] Antoni Viros i Martin, Brian Vaughan, Davis Wertheimer, Joshua Rosenkranz, Mudhakar Srivatsa, Nelson Mimura Gonzalez, Raghu Ganti, Supriyo Chakraborty, Zhuoran Liu, Geeta Chauhan, and Hamid Shojanazeri.  Pytorch compile to speed up inference on llama 2.  PyTorch Blog, 11 2023.
  * [5] Yaniv Leviathan, Matan Kalman, and Yossi Matias.  Fast inference from transformers via speculative decoding.  In Proceedings of the 40th International Conference on Machine Learning, ICML’23. JMLR.org, 2023.
  * [6] Yuhui Li, Fangyun Wei, Chao Zhang, and Hongyang Zhang.  Eagle: Speculative sampling requires rethinking feature uncertainty.  arXiv preprint arXiv:2401.15077, 2024.
  * [7] Xupeng Miao, Gabriele Oliaro, Zhihao Zhang, Xinhao Cheng, Zeyu Wang, Zhengxin Zhang, Rae Ying Yee Wong, Alan Zhu, Lijie Yang, Xiaoxiang Shi, Chunan Shi, Zhuoming Chen, Daiyaan Arfeen, Reyna Abhyankar, and Zhihao Jia.  Specinfer: Accelerating generative large language model serving with tree-based speculative inference and verification.  arXiv preprint arXiv:2305.09781, 2023.
  * [8] Sen Yang, Shujian Huang, Xinyu Dai, and Jiajun Chen.  Multi-candidate speculative decoding.  arXiv preprint arXiv:2401.06706, 2024.
  * [9] Aonan Zhang, Chong Wang, Yi Wang, Xuanyu Zhang, and Yunfei Cheng.  Recurrent drafter for fast speculative decoding in large language models.  arXiv preprint arXiv:2403.09919, 2024.


Assume $u\in C([0,T],L^{1}_{\mathrm{loc}}({\mathbb{R}}^{d}))\cap
L^{\infty}([0,T]\times{\mathbb{R}}^{d})$ is continuous almost everywhere and
$u(T,\cdot)=u_{T}\in L^{\infty}({\mathbb{R}}^{d})$. Then $u$ is the duality
solution of (4.2) if and only if there exist $r^{\delta},r_{\delta}\in
L^{1}_{\mathrm{loc}}([0,T]\times{\mathbb{R}}^{d}))$ such that $\lim_{\delta\to
0}r^{\delta}=\lim_{\delta\to 0}r_{\delta}=0$ in $L^{1}_{\mathrm{loc}}$, and
the $\sup$\- and $\inf$-convolutions

$u^{\delta}(t,x):=\operatorname*{ess\,sup}_{y\in{\mathbb{R}}^{d}}\left\\{u(t,y)-\frac{1}{2\delta}|x-y|^{2}\right\\}$

and

$u_{\delta}(t,x):=\operatorname*{ess\,inf}_{y\in{\mathbb{R}}^{d}}\left\\{u(t,y)+\frac{1}{2\delta}|x-y|^{2}\right\\}$

satisfy in the sense of distributions on $[0,T]\times{\mathbb{R}}^{d}$ the
inequalities

$\frac{\partial u^{\delta}}{\partial t}+b(t,x)\cdot\nabla u^{\delta}\leq
r^{\delta}(t,x)\quad\text{and}\quad\frac{\partial u_{\delta}}{\partial
t}+b(t,x)\cdot\nabla u_{\delta}\geq-r_{\delta}(t,x).$

###### Proof.

Assume first that the $\sup$\- and $\inf$-convolutions have the stated
properties. For standard mollifiers $(\rho_{\eta})_{\eta>0}$ on
${\mathbb{R}}$, define
$u^{\delta}_{\eta}(t,x)=(u^{\delta}(\cdot,x)*_{t}\rho_{\eta})(t)$ and
$u_{\delta,\eta}(t,x)=(u_{\delta}(\cdot,x)*_{t}\rho_{\eta})(t)$. Then, by
Lemma 4.3, $u^{\delta}_{\eta}$ and $u_{\delta,\eta}$ are Lipschitz continuous
on $[0,T]\times{\mathbb{R}}^{d}$, and satisfy a.e. in
$[0,T]\times{\mathbb{R}}^{d}$

$\frac{\partial u^{\delta}_{\eta}}{\partial t}+b(t,x)\cdot\nabla
u^{\delta}_{\eta}\leq r^{\delta}_{\eta}(t,x)\quad\text{and}\quad\frac{\partial
u_{\delta,\eta}}{\partial t}+b(t,x)\cdot\nabla u_{\delta,\eta}\geq-
r_{\delta,\eta}(t,x),$

where

$r^{\delta}_{\eta}(t,x)=(r^{\delta}(\cdot,x)*_{t}\rho_{\eta})(t)+\int_{\mathbb{R}}(b(t,x)-b(s,x))\cdot\nabla
u^{\delta}(s,x)\rho_{\eta}(s-t)ds$

and

$r_{\delta,\eta}(t,x)=(r_{\delta}(\cdot,x)*_{t}\rho_{\eta})(t)+\int_{\mathbb{R}}(b(t,x)-b(s,x))\cdot\nabla
u_{\delta}(s,x)\rho_{\eta}(s-t)ds.$

The (local) boundedness of $b$, $\nabla u_{\delta}$, and $\nabla u^{\delta}$
then allow us to invoke the dominated convergence theorem to say that, for
fixed $\delta$, $\lim_{\eta\to 0}r^{\delta}_{\eta}=r^{\delta}$ and
$\lim_{\eta\to 0}r_{\delta,\eta}=r_{\delta}$ in $L^{1}_{\mathrm{loc}}$.

Now let $f_{0}\in C_{c}({\mathbb{R}}^{d})$ be nonnegative and let $f$ be the
“good” solution of (4.1). In view of the nonnegativity of $J$, $f$ given by
(4.3) is nonnegative on $[0,T]\times{\mathbb{R}}^{d}$, and the bounds for the
backward flow in Lemma 2.2 imply that $f$ has compact support in
$[0,T]\times{\mathbb{R}}^{d}$. By Theorem 4.1, $f$ is a distributional
solution, and therefore

$\displaystyle\int_{{\mathbb{R}}^{d}}f(T,x)u^{\delta}_{\eta}(T,x)dx-\int_{{\mathbb{R}}^{d}}f_{0}(x)u^{\delta}_{\eta}(0,x)dx$
$\displaystyle=\int_{0}^{T}\int_{{\mathbb{R}}^{d}}f(t,x)\left[\partial_{t}u^{\delta,\eta}(t,x)+b(t,x)\cdot\nabla
u^{\delta,\eta}(t,x)\right]dxdt$
$\displaystyle\leq\int_{0}^{T}\int_{{\mathbb{R}}^{d}}f(t,x)r^{\delta}_{\eta}(t,x)dxdt.$

Sending first $\eta\to 0$ and then $\delta\to 0$, using Lemma 4.3 and the
dominated convergence theorem, we conclude that

$\int_{{\mathbb{R}}^{d}}f(T,x)u_{T}(x)dx\leq\int_{{\mathbb{R}}^{d}}f_{0}(x)u(0,x)dx.$

Arguing similarly with $u_{\delta,\eta}$ as a test function, we achieve the
opposite inequality. By linearity, the duality identity holds for any
$f_{0}\in L^{\infty}$ with bounded support, and we conclude that $u$ is the
unique duality solution.

Assume now conversely that $u$ is the duality solution. Let
$(b^{\varepsilon})_{\varepsilon>0}$ be as in (2.3), let $u^{\varepsilon}$ be
the corresponding solution, and define

$u^{\varepsilon,\delta}(t,x):=\sup_{y\in{\mathbb{R}}^{d}}\left\\{u^{\varepsilon}(t,y)-\frac{1}{2\delta}|x-y|^{2}\right\\}$

and

$u^{\varepsilon}_{\delta}(t,x):=\inf_{y\in{\mathbb{R}}^{d}}\left\\{u^{\varepsilon}(t,y)+\frac{1}{2\delta}|x-y|^{2}\right\\}.$

By Lemma 4.3, for fixed $\delta>0$, $u^{\varepsilon,\delta}$ and
$u^{\varepsilon}_{\delta}$ are Lipschitz continuous in the space variable,
uniformly over $[0,T]\times{\mathbb{R}}^{d}$ and $\varepsilon>0$. Moreover,
the $\sup$ and $\inf$ are actually a $\max$ and $\min$, and may be restricted
to

$|y-x|\leq(\max u_{0}-\min u_{0})^{1/2}\delta^{1/2}$

(note that we have used the maximum principle for the transport equation to
control the maximum and minimum of $u^{\varepsilon}$ and $u_{\varepsilon}$).
We may alternatively restrict the $y$ for which the maximum in the definition
of $u^{\varepsilon,\delta}(t,x)$ is attained to satisfy

(4.11) $|y-x|\leq
2(u^{\varepsilon,2\delta}(t,x)-u^{\varepsilon,\delta}(t,x))^{1/2}\delta^{1/2},$

and the minimum in the definition of $u^{\varepsilon}_{\delta}$ is attained by
$y$ satisfying

(4.12) $|y-x|\leq
2(u^{\varepsilon}_{\delta}(t,x)-u^{\varepsilon}_{2\delta}(t,x))^{1/2}\delta^{1/2}.$

Standard properties of envelopes then give the identities, for any
$(t,x)\in[0,T]\times{\mathbb{R}}^{d}$,

$\frac{\partial u^{\varepsilon,\delta}}{\partial t}(t,x)=\frac{\partial
u^{\varepsilon}}{\partial t}(t,y)\quad\text{and}\quad\nabla
u^{\varepsilon,\delta}(t,x)=\nabla u^{\varepsilon}(t,y)=\frac{y-x}{\delta}$

for some $y$ satisfying (4.11). Therefore

$\partial_{t}u^{\varepsilon,\delta}(t,x)=-b^{\varepsilon}(t,y)\cdot\nabla
u^{\varepsilon,\delta}(t,x),$

from which we deduce that $u^{\varepsilon,\delta}$ is uniformly Lipschitz
continuous in the time variable over $[0,T]\times B_{R}$ for any $R>0$,
independently of $\varepsilon$. Further developing the equality gives

(4.13) $\begin{split}\frac{\partial u^{\varepsilon,\delta}}{\partial
t}(t,x)+b^{\varepsilon}(t,x)\cdot\nabla
u^{\varepsilon,\delta}(t,x)&=\frac{\partial u^{\varepsilon}}{\partial
t}(t,y)+b^{\varepsilon}(t,x)\cdot\nabla u^{\varepsilon}(t,y)\\\
&=-(b^{\varepsilon}(t,x)-b^{\varepsilon}(t,y))\cdot\frac{x-y}{\delta}\\\ &\leq
C_{0}(t)\frac{|x-y|^{2}}{\delta}\leq
4C_{0}(t)(u^{\varepsilon,2\delta}(t,x)-u^{\varepsilon,\delta}(t,x)).\end{split}$

We similarly have that $u^{\varepsilon}_{\delta}$ is Lipschitz continuous in
the time variable, locally in space, uniformly over $\varepsilon>0$, and

(4.14) $\frac{\partial u^{\varepsilon}_{\delta}}{\partial
t}(t,x)+b^{\varepsilon}(t,x)\cdot\nabla
u^{\varepsilon}_{\delta}(t,x)\geq-4C_{0}(t)(u^{\varepsilon}_{\delta}(t,x)-u^{\varepsilon}_{2\delta}(t,x)).$

We now claim that, as $\varepsilon\to 0$, $u^{\varepsilon,\delta}$ and
$u^{\varepsilon}_{\delta}$ converge pointwise to respectively $u^{\delta}$ and
$u_{\delta}$, and then, by the uniform-in-$\varepsilon$ Lipschitz regularity,
the convergence is locally uniform. To see this, fix $x\in{\mathbb{R}}^{d}$
and $\eta>0$, and let $A\subset{\mathbb{R}}^{d}$ be a set of positive measure
such that

$u^{\delta}(t,x)\leq u(t,y)-\frac{|x-y|^{2}}{2\delta}+\eta.$

We then have, for all $y\in A$,

$u^{\delta,\varepsilon}(t,x)\geq
u^{\varepsilon}(t,y)-\frac{|x-y|^{2}}{2\delta}.$

For at least one such $y$, we then have
$u^{\varepsilon}(t,y)\xrightarrow{\varepsilon\to 0}u(t,y)$, and we thus have

$\limsup_{\varepsilon\to
0}\left(u^{\delta}(t,x)-u^{\delta,\varepsilon}(t,x)\right)\leq\eta.$

It follows that $\limsup_{\varepsilon\to
0}\left(u^{\delta}(t,x)-u^{\delta,\varepsilon}(t,x)\right)\leq 0$ since $\eta$
was arbitrary.

Now, there exists a full measure set $B\subset{\mathbb{R}}^{d}$ such that, for
all $y\in B$,

$u^{\delta}(t,x)\geq
u(t,y)-\frac{|x-y|^{2}}{2\delta}\quad\text{and}\quad\lim_{\varepsilon\to
0}u^{\varepsilon}(t,y)=u(t,y).$

In view of the continuity of $u^{\varepsilon}(t,\cdot)$, there exists a
bounded (independently $\varepsilon$) sequence
$(y_{n})_{n\in{\mathbb{N}}}\subset B$ such that

$\rho_{n}:=u^{\delta,\varepsilon}(t,x)-\left\\{u^{\varepsilon}(t,y_{n})-\frac{|x-y_{n}|^{2}}{2\delta}\right\\}$

satisfies $\lim_{n\to\infty}\rho_{n}=0$. Therefore, for all $n$,

$u^{\delta,\varepsilon}(t,x)-u^{\delta}(t,x)\leq
u^{\varepsilon}(t,y_{n})-u(t,y_{n})+\rho_{n}.$

Sending $\varepsilon\to 0$ gives $\limsup_{\varepsilon\to
0}(u^{\delta,\varepsilon}(t,x)-u^{\delta}(t,x))\leq\rho_{n}$, and the proof of
pointwise convergence is finished upon sending $n\to\infty$. The exact same
argument can be used for the pointwise convergence of
$u^{\varepsilon}_{\delta}$ to $u_{\delta}$.

It then follows that, for fixed $\delta$, as $\varepsilon\to 0$, $\nabla
u^{\varepsilon,\delta}$ and $\nabla u^{\varepsilon}_{\delta}$ converge
weak-$\star$ in $L^{\infty}$ to $\nabla u^{\delta}$ and $\nabla u_{\delta}$
respectively, while $b^{\varepsilon}$ converges in $L^{1}_{\mathrm{loc}}$ to
$b$. We may then take $\varepsilon\to 0$ in (4.13) and (4.14) to obtain the
distributional inequalities

$\frac{\partial u^{\delta}}{\partial t}(t,x)+b(t,x)\cdot\nabla
u^{\delta}(t,x)\leq
4C_{0}(t)(u^{2\delta}(t,x)-u^{\delta}(t,x))=:r^{\delta}(t,x)$

and

$\frac{\partial u_{\delta}}{\partial t}(t,x)+b(t,x)\cdot\nabla
u_{\delta}(t,x)\geq-4C_{0}(t)(u_{\delta}(t,x)-u_{2\delta}(t,x))=:-r_{\delta}(t,x).$

By Lemma 4.3 and the almost-everywhere continuity of $u$, the right-hand sides
of both inequalities converge a.e. to $0$ as $\delta\to 0$, and, by the
uniform boundedness in $\delta$ of $u^{\delta}$ and $u_{\delta}$ and the
dominated convergence theorem, $r^{\delta}$ and $r_{\delta}$ both converge in
$L^{1}_{\mathrm{loc}}$ to $0$ as $\delta\to 0$. ∎

####  4.4.2 The conservative equation: uniqueness of nonnegative solutions

We observe that, in the first implication in the proof of Theorem 4.10, it was
proved that $u$ was a duality solution by proving the duality identity
relative to the a “good” nonnegative solution. However, it was only explicitly
used that $f$ was a distributional solution. Therefore, after having proved
the equivalence in Theorem 4.10, we arrive at the following.

######  Theorem 4.11.

Suppose that $f\in C([0,T],L^{p}_{\mathrm{loc}}({\mathbb{R}}^{d}))$ is a
distributional solution of (4.1) and $f\geq 0$. Then
$f(t,x)=f(0,\phi_{0,t}(x))J_{0,t}(x)$.

###### Proof.

Fix $t>0$ and $v\in C_{c}({\mathbb{R}}^{d})$, and let $u\in
C([0,t],L^{1}_{\mathrm{loc}}({\mathbb{R}}^{d}))\cap
L^{\infty}([0,t]\times{\mathbb{R}}^{d})$ be the duality solution of (4.2) with
terminal data $v$ at time $t$. Then, by Theorem 4.7, $u$ is continuous almost
everywhere in $[0,t]\times{\mathbb{R}}^{d}$. Arguing exactly as in the first
part of Theorem 4.10, using the nonnegativity of $f$, we arrive at the
equality

$\int_{{\mathbb{R}}^{d}}f(t,x)v(x)dx=\int_{{\mathbb{R}}^{d}}f(0,x)u(0,x)dx.$

Since $v$ was arbitrary, it follows from the definition of duality solutions
that $f(t,x)$ must be given by (4.3). ∎

We then have the following corollary about characterizing the good solution
even when $f$ is signed:

######  Corollary 4.3.

A function $f\in C([0,T],L^{p}_{\mathrm{loc}}({\mathbb{R}}^{d}))$ is the good
solution of (4.1) if and only if $f$ and $|f|$ are both solutions in the sense
of distributions.

###### Proof.

That this property is satisfied by the good solution was already pointed out
(Corollary 4.1). Suppose now that $f$ and $|f|$ are both distributional
solutions. It follows that $f_{+}=\frac{1}{2}(f+|f|)$ and
$f_{-}=\frac{1}{2}(|f|-f)$ are distributional solutions, and, since $f_{+}\geq
0$ and $f_{-}\geq 0$, they are both the good solutions. Therefore
$f=f_{+}-f_{-}$ is a good solution by linearity. ∎

####  4.4.3 Uniqueness of regular Lagrangian flows

We can finally establish the uniqueness for the forward flows of the ODE (2.2)

######  Theorem 4.12.

For every $s\in[0,T]$ and almost every $x\in{\mathbb{R}}^{d}$, $\phi_{t,s}(x)$
is the unique absolutely continuous solution of (2.2).

###### Proof.

This is a consequence of Theorem 4.11 and the superposition principle of
Ambrosio [5, Theorem 3.1]. ∎

###  4.5 Some remarks for second order equations

We next investigate the second-order analogues of (4.1) and (4.2). As
mentioned earlier, we are not able to treat the most general case in which
$\sigma$ is a regular function of $x$. This is due to the fact that Lemma 2.5
only gives regularity of the backward stochastic flow in $C^{0,1-\varepsilon}$
for $0<\varepsilon<1$. As a consequence, defining the Jacobian and using it to
analyze the right-inverse of the flow is not possible in general. Our results
in this case are limited to stochastic flows for which the coefficient
$\sigma$ in front of the Wiener process is constant in the space variable. The
generalization to regular but nonconstant $\sigma$ will be the subject of
future work.

####  4.5.1 The expansive stochastic flow with constant noise coefficient

The stochastic analogue of the forward flow (2.2) is

(4.15)
$d_{t}\Phi_{t,s}(x)=b(t,\Phi_{t,s}(x))dt+\sigma(t,\Phi_{t,s}(x))dW_{t},\quad
t\in[s,T],\quad\Phi_{s,s}(x)=x,$

where $\sigma:[0,T]\times{\mathbb{R}}^{d}\to{\mathbb{R}}^{d\times m}$ is some
matrix-valued map. As we shall see, this general setting is out of the reach
at the moment, and we thus assume

(4.16) $\sigma\in L^{2}([0,T],{\mathbb{R}}^{d\times m})$

is constant in the space variable. We then consider the forward stochastic
flow

(4.17) $d\Phi_{t,s}(x)=b(t,\Phi_{t,s}(x))dt+\sigma_{t}dW_{t},\quad
t\in[s,T],\quad\Phi_{s,s}(x)=x.$

Formally defining

$\tilde{\Phi}_{t,s}(x):=\Phi_{t,s}(x)-\underbrace{\int_{s}^{t}\sigma_{r}dW_{r}}_{:=M_{t}-M_{s}}$

leads to the random ODE

(4.18)
$\partial_{t}\tilde{\Phi}_{t,s}(x)=b\left(t,\tilde{\Phi}_{t,s}(x)+M_{t}-M_{s}\right),\quad
t\in[s,T],\quad\tilde{\Phi}_{s,s}(x)=x.$

We now invoke the theory of the previous subsections to obtain the following:

######  Theorem 4.13.

For every $s\in[0,T)$, with probability one, there exists a unique
$\Phi_{\cdot,s}\in C([s,T],L^{p}_{\mathrm{loc}}({\mathbb{R}}^{d}))\cap
L^{p}_{\mathrm{loc}}({\mathbb{R}}^{d},C([s,T]))$ such that, for a.e.
$x\in{\mathbb{R}}^{d}$,

$\Phi_{t,s}(x)=x+\int_{s}^{t}b(r,\Phi_{r,s}(x))dr+\int_{s}^{t}\sigma_{r}dW_{r}.$

If $(b^{\varepsilon})_{\varepsilon>0}$ are as in (2.3) and
$\Phi^{\varepsilon}$ is the unique stochastic flow solving (4.17) with drift
$b^{\varepsilon}$, then, with probability one, as $\varepsilon\to 0$,
$\Phi^{\varepsilon}$ converges in
$C([s,T],L^{p}_{\mathrm{loc}}({\mathbb{R}}^{d}))$ and in
$L^{p}_{\mathrm{loc}}({\mathbb{R}}^{d},C([s,T]))$ to $\Phi$.

###### Proof.

This follows upon applying the results of Theorems 4.8 and 4.12 to the random
ODE (4.18). ∎

####  4.5.2 A priori estimates for the second-order nonconservative equation

We next relate the forward stochastic flow from the previous subsection to the
terminal value problem for a certain second-order, nonconservative equation.
This will be done with the use of a priori $L^{p}$ and $BV$ estimates, which
lead to useful compactness results, just as for the first order case.

We begin with the more general problem

(4.19) $-\partial_{t}u-\operatorname{tr}[a(t,x)\nabla^{2}u]+b(t,x)\cdot\nabla
u=0\quad\text{in }(0,T)\times{\mathbb{R}}^{d},\quad u(T,\cdot)=u_{T},$

where

(4.20) $a(t,x)=\frac{1}{2}\sigma(t,x)\sigma(t,x)^{T},\quad\sigma\in
L^{2}([0,T],C^{1,1}({\mathbb{R}}^{d},{\mathbb{R}}^{d\times m}));$

notice that, although we allow $\sigma$ to be nonconstant here, we require
more regularity for $\sigma$ than in Section 3.

######  Lemma 4.4.

There exists $C\in L^{1}_{+}([0,T])$ depending only the $C^{1,1}$ norm of
$\sigma$ such that, if $u$ is a smooth solution of

$-\partial_{t}u-\operatorname{tr}[a(t,x)\nabla^{2}u]=0\quad\text{in
}(0,T)\times{\mathbb{R}}^{d},\quad u(T,\cdot)=u_{T},$

then

$\left\|u(t,\cdot)\right\|_{BV({\mathbb{R}}^{d})}\leq\exp\left(\int_{t}^{T}C(s)ds\right)\left\|u_{T}\right\|_{BV({\mathbb{R}}^{d})}.$

###### Proof.

For $(t,x,z)\in[0,T]\times{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}$, set
$w(t,x,z)=\nabla u(t,x)\cdot z$. Then $w$ solves the parabolic PDE

$\frac{\partial w}{\partial
t}-\operatorname{tr}[A(t,x,z)\nabla^{2}_{(x,z)}w]=0\quad\text{in
}(0,T)\times{\mathbb{R}}^{2d},$

where

$A(t,x,z)=\frac{1}{2}\begin{pmatrix}\sigma(t,x)\\\
z\cdot\nabla\sigma(t,x)\end{pmatrix}\begin{pmatrix}\sigma(t,x)^{T}&z\cdot\nabla\sigma(t,x)^{T}\end{pmatrix}.$

After a routine regularization argument, using the convexity of $w\mapsto|w|$,

(4.21) $\frac{\partial|w|}{\partial
t}-\operatorname{tr}[A(t,x,z)\nabla^{2}_{(x,z)}|w|]\leq 0\quad\text{in
}(0,T)\times{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}.$

For some $m>d+1$, let $\phi\in C^{\infty}_{+}([0,\infty))$ be such that, for
some universal $C>0$,

(4.22) $\phi(r)=\frac{1}{r^{m}}\quad\text{for }r\geq 1\quad\text{and}\quad
r|\phi^{\prime}(r)|+r^{2}|\phi^{\prime\prime}(r)|\leq C\phi(r)\quad\text{for
all }r\geq 0.$

We multiply (4.21) by $\phi(|z|)$ and integrate in
$(x,z)\in{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}$. Then (4.20) and (4.22) imply
that for some $C\in L^{1}_{+}([0,T])$,

$-\frac{d}{dt}\iint_{{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}}|w(t,x,z)|\phi(|z|)dxdz\leq
C(t)\iint_{{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}}|w(t,x,z)|\phi(|z|)dxdz.$

The proof is then finished by Grönwall’s lemma and the fact that

$\iint_{{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}}|w(t,x,z)|\phi(z)dxdz=c_{0}\int_{{\mathbb{R}}^{d}}|\nabla
u(t,x)|dx,$

where $c_{0}:=\int_{{\mathbb{R}}^{d}}|\nu\cdot z|\phi(|z|)dz$ is finite and
independent of $|\nu|=1$. ∎

We have already proved an exponential propagation of the $BV$ bounds when
$a=0$ in Lemma 4.1. It is a classical fact for evolution PDEs that, upon using
a splitting scheme, that these estimates can be combined, and we immediately
have the following:

######  Lemma 4.5.

There exists a constant $C\in L^{1}_{+}([0,T])$ depending only on the
constants in (2.1) and (4.20) such that, if $u$ is a smooth solution of
(4.19), then

$\left\|u(t,\cdot)\right\|_{L^{p}}\leq\exp\left(\int_{0}^{t}C(s)ds\right)\left\|u_{T}\right\|_{L^{p}}\quad\text{and}\quad\left\|u(t,\cdot)\right\|_{BV}\leq\exp\left(\int_{0}^{t}C(s)ds\right)\left\|u_{T}\right\|_{BV}.$

Just as in the first-order case, it is not possible to define
$L^{p}$-distributional solutions of (4.19), and the utility of Lemma 4.5 is
that it allows to obtain strongly convergent subsequences in
$C([0,T],L^{p}({\mathbb{R}}^{d}))$ after regularizing the velocity field $b$.

The main question is whether such limiting solutions are unique. This
uniqueness was achieved in the first-order case through duality with the
conservative equation, and the solution was further characterized with a
formula involving the forward flow. In the second-order case, we are
constrained to work with constant noise coefficients:

(4.23) $-\partial_{t}u-\operatorname{tr}[a(t)\nabla^{2}u]+b(t,x)\cdot\nabla
u=0\quad\text{in }(0,T)\times{\mathbb{R}}^{d},\quad u(T,\cdot)=u_{T},$

where $a=\frac{1}{2}\sigma\sigma^{T}$ as before.

######  Theorem 4.14.

For $1<p<\infty$ and $t\in[0,T]$, the map

$C_{c}({\mathbb{R}}^{d})\ni u_{T}\mapsto{\mathbb{E}}[u_{T}\circ\Phi_{T,t}]$

extends to a continuous, linear, order-preserving map on
$L^{p}({\mathbb{R}}^{d})$, and the function

(4.24)
$u(t,x):={\mathbb{E}}[u_{T}(\Phi_{T,t}(x))]\quad(t,x)\in[0,T]\times{\mathbb{R}}^{d}$

belongs to $C([0,T],L^{p}({\mathbb{R}}^{d}))$, and, if $u_{T}\in
BV({\mathbb{R}}^{d})$, then $u\in L^{\infty}([0,T],BV({\mathbb{R}}^{d}))$.

If $(b^{\varepsilon})_{\varepsilon>0}$ is as in (2.3) and $u^{\varepsilon}$ is
the corresponding solution of (4.23), then, as $\varepsilon\to 0$,
$u^{\varepsilon}$ converges strongly to $u$ in
$C([0,T],L^{p}({\mathbb{R}}^{d}))$.

###### Proof.

Assume that $u_{T}\in C^{2}({\mathbb{R}}^{d})\cap C_{c}({\mathbb{R}}^{d})$.
For $b^{\varepsilon}$ and $u^{\varepsilon}$ as in the statement of the
theorem, we have the standard representation formula
$u^{\varepsilon}(t,x)={\mathbb{E}}[u_{T}(\Phi^{\varepsilon}_{T,t}(x))]$, where
$\Phi^{\varepsilon}$ corresponds to the flow (4.17) with drift
$b^{\varepsilon}$. By Theorem 4.13, for any $t\in[0,T]$, with probability one,
$u_{T}\circ\Phi^{\varepsilon}_{T,t}\to u_{T}\circ\Phi_{T,t}$ a.e. in
${\mathbb{R}}^{d}$. On the other hand, by Lemma 4.5,
$(u^{\varepsilon})_{\varepsilon>0}$ is precompact in
$C([0,T],L^{p}({\mathbb{R}}^{d}))$, and therefore the full sequence converges
to $u$ given by (4.24). The $L^{p}$-bounds and the extension to $u_{T}\in
L^{p}({\mathbb{R}}^{d})$ now follow from the $L^{p}$ a priori estimates in
Lemma 4.5. ∎

####  4.5.3 Representation formula for the Fokker-Planck equation

We turn next to the Fokker-Planck equation

(4.25)
$\partial_{t}f-\nabla^{2}\cdot(a(t,x)f)+\operatorname{div}(b(t,x)f)=0\quad\text{in
}(0,T)\times{\mathbb{R}}^{d},\quad f(0,\cdot)=f_{0},$

where once again $a=\frac{1}{2}\sigma\sigma^{T}$ with $\sigma$ as in (4.20).

The existence of solutions in $C([0,T],L^{p}({\mathbb{R}}^{d}))$ is
straightforward; we include the proof for convenience.

######  Theorem 4.15.

For any $f_{0}\in L^{p}({\mathbb{R}}^{d})$, $1\leq p\leq\infty$, there exists
a distributional solution $f\in C([0,T],L^{p}_{\mathrm{w}}({\mathbb{R}}^{d}))$
if $1\leq p<\infty$, or $f\in L^{\infty}$ if $p=\infty$. Moreover, there
exists $C\in L^{1}_{+}([0,T])$ depending only on $p$, $C_{0}(t)$ from (2.1)
and the $L^{2}([0,T],C^{1,1}({\mathbb{R}}^{d}))$ norm of $a$888In fact, only
an upper bound for $\nabla^{2}\cdot a=\partial_{ij}a_{ij}$ is needed. such
that

$\left\|f(t,\cdot)\right\|_{L^{p}}\leq\exp\left(\int_{0}^{t}C(s)ds\right)\left\|f\right\|_{L^{p}}.$

###### Proof.

We do this with the use of a priori estimates, assuming all the data is
smooth. The computations may be made rigorous by regularizing $b$, adding a
small ellipticity to $a$, and extracting weakly convergent subsequences.

We then compute

$\partial_{t}|f|^{p}-\nabla^{2}\cdot(a(t,x)|f|^{p})+\operatorname{div}(b(t,x)|f|^{p})\leq(p-1)\left(\nabla^{2}\cdot
a(t,x)-\operatorname{div}b(t,x)\right)|f|^{p},$

and so $\partial_{t}\int|f(t,\cdot)|^{p}\leq C(t)\int|f(t,\cdot)^{p}$ for some
$C$ as in the statement of the Theorem. The result now follows from Grönwall’s
lemma. ∎

We now explore the possibility of obtaining a formula for the solution,
similar to (4.3) for the first order equation (4.1). To do so, it is
convenient to reverse time and consider, for fixed $t\in(0,T]$, the equation
satisfied by $g^{(t)}(s,x):=f(t-s,x)$:

$-\partial_{s}g^{(t)}-\nabla^{2}\cdot(a(t-s,x)g^{(t)})+\operatorname{div}(b(t-s,x)g^{(t)})=0\quad\text{in
}(0,t)\times{\mathbb{R}}^{d},\quad g^{(t)}(t,\cdot)=f_{0}.$

For $(s,x,\xi)\in[0,t]\times{\mathbb{R}}^{d}\times{\mathbb{R}}$, define
$G^{(t)}(s,x,\xi)=g^{(t)}(s,x)\xi$. Then

(4.26)
$\begin{dcases}-\partial_{s}G^{(t)}-\operatorname{tr}[A^{(t)}(s,x,\xi)\nabla^{2}_{x,\xi}G^{(t)}]-B^{(t)}(s,x)\cdot\nabla
G^{(t)}-C^{(t)}(s,x)\xi\partial_{\xi}G^{(t)}=0&\text{in
}(0,t)\times{\mathbb{R}}^{d+1},\\\ G^{(t)}(t,x,\xi)=f_{0}(x)\xi,\end{dcases}$

where

(4.27)
$\left\\{\begin{split}A^{(t)}(s,x,\xi)&=\frac{1}{2}\Sigma^{(t)}(s,x,\xi)\Sigma^{(t)}(s,x,\xi)^{T},\quad\Sigma^{(t)}(s,x,\xi)=\begin{pmatrix}\sigma\\\
\xi\operatorname{div}\sigma\end{pmatrix},\\\
B^{(t)}(s,x)&=-b+(\sigma\cdot\nabla)\sigma^{T},\quad\text{and}\\\
C^{(t)}(s,x)&=-\operatorname{div}\left(b-\operatorname{div}a\right)\\\
&=-\operatorname{div}b+\operatorname{tr}[(\sigma\cdot\nabla)(\nabla\cdot\sigma)]+\frac{1}{2}|\operatorname{div}\sigma|^{2}+\frac{1}{2}\operatorname{tr}[\nabla\sigma\nabla\sigma^{T}];\end{split}\right.$

for brevity, we have suppressed the arguments for $a$, $\sigma$, and $b$,
which are all $(t-s,x)$.

For an $m$-dimensional Wiener process $W$ on $[0,t]$ and a fixed $s\in[0,t]$,
we are led to consider the SDE, for $r\in[s,t]$,

(4.28) $\begin{dcases}d_{r}\begin{pmatrix}\Phi^{(t)}_{r,s}(x,\xi)\\\
\Xi^{(t)}_{r,s}(x,\xi)\end{pmatrix}=\begin{pmatrix}B^{(t)}(r,\Phi^{(t)}_{r,s}(x,\xi))\\\
C^{(t)}(r,\Phi^{(t)}_{r,s}(x,\xi))\Xi^{(t)}_{r,s}(x,\xi)\end{pmatrix}dr+\Sigma^{(t)}(r,\Phi^{(t)}_{r,s}(x,\xi),\Xi^{(t)}_{r,s}(x,\xi))dW_{r},\\\
\begin{pmatrix}\Phi^{(t)}_{s,s}(x,\xi)\\\
\Xi^{(t)}_{s,s}(x,\xi)\end{pmatrix}=\begin{pmatrix}x\\\
\xi\end{pmatrix}.\end{dcases}$

Itô’s formula, (4.26), and (4.28) then yield that, for any
$(s,x,\xi)\in[0,t)\times{\mathbb{R}}^{d}\times{\mathbb{R}}$,

$r\mapsto G^{(t)}(r,\Phi^{(t)}_{r,s}(x,\xi),\Xi^{(t)}_{r,s}(x,\xi))$

is a martingale on $[s,t]$ with respect to the filtration
$(\mathcal{F}_{r})_{r\in[0,t]}$ generated by the Wiener process $W$, and so,
for all $r\in[s,t]$,

(4.29)
${\mathbb{E}}\left[G^{(t)}(r,\Phi^{(t)}_{r,s}(x,\xi),\Xi^{(t)}_{r,s}(x,\xi))\mid\mathcal{F}_{s}\right]=G^{(t)}(s,x,\xi).$

Observe that $\Phi^{(t)}_{r,s}$ is independent of $\xi$, while
$\Xi^{(t)}_{r,s}$ can be written as
$\Xi^{(t)}_{r,s}(x,\xi)=J^{(t)}_{r,s}(x)\xi$ for some scalar quantity
$J^{(t)}_{r,s}(x)$, and so (4.28) reduces to the two SDEs

(4.30)
$\begin{dcases}d_{r}\Phi^{(t)}_{r,s}(x)=-\left[b(t-r,\Phi^{(t)}_{r,s}(x))-(\sigma\cdot\nabla)\sigma^{T}(t-r,\Phi^{(t)}_{r,s}(x))\right]dt+\sigma(t-r,\Phi^{(t)}_{r,s}(x))dW_{r},\quad
r\in[s,t],\\\ \Phi^{(t)}_{s,s}(x)=x\end{dcases}$

and

(4.31)
$\left\\{\begin{split}&d_{r}J^{(t)}_{r,s}(x)=\left[-\operatorname{div}b+\operatorname{tr}[(\sigma\cdot\nabla)(\nabla\cdot\sigma)]+\frac{1}{2}|\operatorname{div}\sigma|^{2}+\frac{1}{2}\operatorname{tr}[\nabla\sigma\nabla\sigma^{T}]\right](t-r,\Phi^{(t)}_{r,s}(x))J^{(t)}_{r,s}(x)dr\\\
&\qquad+\operatorname{div}\sigma(t-r,\Phi^{(t)}_{r,s}(x))J^{(t)}_{r,s}(x)dW_{r},\quad
r\in[s,t],\\\ &J^{(t)}_{s,s}(x)=1.\end{split}\right.$

Standard but tedious computations involving Itô’s formula reveal that
$J^{(t)}_{r,s}(x)=\det\nabla_{x}\Phi^{(t)}_{r,s}(x)$.

Taking $r=t$ and $\xi=1$ in (4.29), we thus arrive at

${\mathbb{E}}\left[f_{0}(\Phi^{(t)}_{t,s}(x))J^{(t)}_{t,s}(x)\mid\mathcal{F}_{s}\right]=g(s,x),$

and so, because $g(0,x)=f(t,x)$, we obtain the representation for solutions of
(4.25):

(4.32)
$f(t,x)={\mathbb{E}}\left[f_{0}(\Phi^{(t)}_{t,0}(x))J^{(t)}_{t,0}(x)\right].$

Let us note that $\Phi^{(t)}_{t,0}$ has the same law as $(\Phi_{t,0})^{-1}$,
where $\Phi_{t,s}$ is the stochastic flow from (4.15). We can see this by
duality with the nonconservative equation. Indeed, if $u$ is the solution of
(4.19) with $u(t,\cdot)=g$ for some given $g$, then

$\int f_{0}(x)u(0,x)dx=\int f(t,x)g(x)dx.$

On the other hand, by (4.24) and (4.32),

$\int f_{0}(x)u(0,x)dx={\mathbb{E}}\int f_{0}(x)g(\Phi_{t,0}(x))dx$

and

$\int f(t,x)g(x)dx={\mathbb{E}}\int
f_{0}(\Phi^{(t)}_{t,0}(x))g(x)J^{(t)}_{t,0}(x)dx,$

so, using the change of variables formula and the fact that $f_{0}$ is
arbitrary, we have
${\mathbb{E}}[g(\Phi_{t,0}(x))]={\mathbb{E}}[g([\Phi^{(t)}_{t,0}]^{-1}(x))]$
for all $g:{\mathbb{R}}^{d}\to{\mathbb{R}}$ and $x\in{\mathbb{R}}^{d}$.

We now note that the SDE (4.30) falls under the assumptions of Lemma 2.5, and
therefore, for every $0\leq s<t\leq T$, there exists a unique solution
$\Phi^{(t)}_{\cdot,s}$ with the properties laid out by that result. However,
the main difficulty is that we do not know whether $\Phi^{(t)}_{t,0}$ is
Lipschitz continuous on ${\mathbb{R}}^{d}$ (see Remark 2.8). This prevents us
from bounding $J^{(t)}_{t,0}$ uniformly in $L^{\infty}$ and passing to weak
distributional limits. This is a major obstacle in using the formula (4.32) to
identify the unique limiting distributional solution of (4.25), as we did for
the first order equation (4.1).

The exception is when $\sigma$ is independent of $x$. In that case, (4.30) and
(4.31) become

(4.33)
$d_{r}\Phi^{(t)}_{r,s}(x)=-b(t-r,\Phi^{(t)}_{r,s}(x))dr+\sigma(t-r)dW_{r},\quad
r\in[s,t],\quad\Phi^{(t)}_{s,s}(x)=x$

and

(4.34)
$\partial_{r}J^{(t)}_{r,s}(x)=-\operatorname{div}b(t-r,\Phi^{(t)}_{r,s}(x))J^{(t)}_{r,s}(x),\quad
r\in[s,t],\quad J^{(t)}_{s,s}(x)=1.$

The SDE (4.34) is in fact an ODE with random coefficients. In particular,
$J^{(t)}_{\cdot,s}$ has a deterministic bound.

We then characterize uniquely the limiting distributional solution of

(4.35)
$\partial_{t}f-\nabla^{2}\cdot(a(t)f)+\operatorname{div}(b(t,x)f)=0\quad\text{in
}(0,T)\times{\mathbb{R}}^{d},\quad f(0,\cdot)=f_{0}.$

######  Theorem 4.16.

For $1\leq p<\infty$, the formula (4.32), where $\Phi^{(t)}_{\cdot,s}$ and
$J^{(t)}_{\cdot,s}$ are specified by respectively (4.33) and (4.34), extends
continuously to any $f_{0}\in L^{p}({\mathbb{R}}^{d})$. If $f_{0}\in
L^{p}({\mathbb{R}}^{d})$ and $(b^{\varepsilon})_{\varepsilon>0}$ are as in
(2.3) and $f^{\varepsilon}$ is the corresponding solution of (4.35), then, as
$\varepsilon\to 0$, $f^{\varepsilon}$ converges weakly in
$C([0,T],L^{p}_{\mathrm{w}}({\mathbb{R}}^{d}))$ to $f$. If $f_{0}\geq 0$, then
there exists a unique nonnegative distributional solution of (4.35), which is
given by (4.32).

###### Proof.

Let $(b^{\varepsilon})_{\varepsilon>0}$ and $f^{\varepsilon}$ be as in the
statement of the theorem, and assume $f_{0}\in C^{2}_{c}({\mathbb{R}}^{d})$.
Let $u^{\varepsilon}$ be the solution of (4.23) with velocity
$b^{\varepsilon}$ and with terminal data $u(t,\cdot)=g\in
C^{2}_{c}({\mathbb{R}}^{d})$ for some fixed $t\in[0,T]$. Then integration by
parts yields

$\int f^{\varepsilon}(t,x)g(x)dx=\int f_{0}(x)u^{\varepsilon}(0,x)dx.$

By Theorem 4.14, as $\varepsilon\to 0$, $u^{\varepsilon}$ converges strongly
in $L^{p^{\prime}}({\mathbb{R}}^{d})$ to the function $u$ defined uniquely by
$u(s,x)=g(\Phi_{t,s}(x))$. Therefore, any
$C([0,T],L^{p}_{\mathrm{w}}({\mathbb{R}}^{d}))$-weak limit $f$ of
$f^{\varepsilon}$ as $\varepsilon\to 0$ must satisfy

$\int f(t,x)g(x)dx=\int f_{0}(x)u(0,x)dx,$

and it follows that there is a unique such limiting function $f$.

On the other hand, for $\varepsilon>0$,

$f^{\varepsilon}(t,x)={\mathbb{E}}\left[f_{0}\left(\Phi^{(t),\varepsilon}_{t,0}(x)\right)J^{(t),\varepsilon}_{t,0}(x)\right],$

where $\Phi^{(t),\varepsilon}_{\cdot,s}$ and $J^{(t),\varepsilon}_{\cdot,s}$
are as in respectively (4.33) and (4.34) with $b$ replaced everywhere by
$b^{\varepsilon}$. For fixed $t\in[0,T]$, uniformly in $\varepsilon$,
$\Phi^{(t),\varepsilon}_{t,0}$ is Lipschitz continuous on ${\mathbb{R}}^{d}$,
and so $J^{(t),\varepsilon}_{t,0}=\det\nabla_{x}\Phi^{(t),\varepsilon}_{t,0}$
is bounded in $L^{\infty}$. By exactly the same arguments as in Lemma 2.3 and
Theorem 4.1, we see that, as $\varepsilon\to 0$,
${\mathbb{E}}f_{0}\left(\Phi^{(t),\varepsilon}_{t,0}\right)J^{(t),\varepsilon}_{t,0}$
converges weakly in $L^{p}$ to
${\mathbb{E}}f_{0}\left(\Phi^{(t)}_{t,0}\right)J^{(t)}_{t,0}$. It follows that
$f$ must be given by (4.32). The fact that the formula extends to arbitrary
$f_{0}\in L^{p}({\mathbb{R}}^{d})$ now follows from the a priori $L^{p}$
bounds in Theorem 4.15.

The uniqueness of nonnegative distributional solutions is then a consequence
of the uniqueness of the forward flow established in Theorem 4.13, as well as
the generalization of superposition to second-order Fokker-Planck equations
(see Figalli [38, Lemma 2.3]). ∎

## References

  * [1] Alberti, G., Bianchini, S., and Crippa, G. Two-dimensional transport equation with Hamiltonian vector fields.  In Hyperbolic problems: theory, numerics and applications, vol. 67, Part 2 of Proc. Sympos. Appl. Math. Amer. Math. Soc., Providence, RI, 2009, pp. 337–346.
  * [2] Alberti, G., Bianchini, S., and Crippa, G. Structure of level sets and Sard-type properties of Lipschitz maps.  Ann. Sc. Norm. Super. Pisa Cl. Sci. (5) 12, 4 (2013), 863–902.
  * [3] Alberti, G., Bianchini, S., and Crippa, G. A uniqueness result for the continuity equation in two dimensions.  J. Eur. Math. Soc. (JEMS) 16, 2 (2014), 201–234.
  * [4] Ambrosio, L. Transport equation and Cauchy problem for $BV$ vector fields.  Invent. Math. 158, 2 (2004), 227–260.
  * [5] Ambrosio, L. Transport equation and Cauchy problem for non-smooth vector fields.  In Calculus of variations and nonlinear partial differential equations, vol. 1927 of Lecture Notes in Math. Springer, Berlin, 2008, pp. 1–41.
  * [6] Ambrosio, L. Well posedness of ODE’s and continuity equations with nonsmooth vector fields, and applications.  Rev. Mat. Complut. 30, 3 (2017), 427–450.
  * [7] Ambrosio, L., Bouchut, F., and De Lellis, C. Well-posedness for a class of hyperbolic systems of conservation laws in several space dimensions.  Comm. Partial Differential Equations 29, 9-10 (2004), 1635–1651.
  * [8] Ambrosio, L., Colombo, M., De Philippis, G., and Figalli, A. Existence of Eulerian solutions to the semigeostrophic equations in physical space: the 2-dimensional periodic case.  Comm. Partial Differential Equations 37, 12 (2012), 2209–2227.
  * [9] Ambrosio, L., Colombo, M., De Philippis, G., and Figalli, A. A global existence result for the semigeostrophic equations in three dimensional convex domains.  Discrete Contin. Dyn. Syst. 34, 4 (2014), 1251–1268.
  * [10] Ambrosio, L., Colombo, M., and Figalli, A. Existence and uniqueness of maximal regular flows for non-smooth vector fields.  Arch. Ration. Mech. Anal. 218, 2 (2015), 1043–1081.
  * [11] Ambrosio, L., and Crippa, G. Continuity equations and ODE flows with non-smooth velocity.  Proc. Roy. Soc. Edinburgh Sect. A 144, 6 (2014), 1191–1244.
  * [12] Ambrosio, L., Crippa, G., and Maniglia, S. Traces and fine properties of a $BD$ class of vector fields and applications.  Ann. Fac. Sci. Toulouse Math. (6) 14, 4 (2005), 527–561.
  * [13] Ambrosio, L., and De Lellis, C. Existence of solutions for a class of hyperbolic systems of conservation laws in several space dimensions.  Int. Math. Res. Not., 41 (2003), 2205–2220.
  * [14] Ambrosio, L., De Lellis, C., and Malý, J. On the chain rule for the divergence of BV-like vector fields: applications, partial results, open problems.  In Perspectives in nonlinear partial differential equations, vol. 446 of Contemp. Math. Amer. Math. Soc., Providence, RI, 2007, pp. 31–67.
  * [15] Ambrosio, L., Lecumberry, M., and Maniglia, S. Lipschitz regularity and approximate differentiability of the DiPerna-Lions flow.  Rend. Sem. Mat. Univ. Padova 114 (2005), 29–50.
  * [16] Ambrosio, L., and Trevisan, D. Lecture notes on the DiPerna-Lions theory in abstract measure spaces.  Ann. Fac. Sci. Toulouse Math. (6) 26, 4 (2017), 729–766.
  * [17] Bianchini, S., and Bonicatto, P. A uniqueness result for the decomposition of vector fields in $\mathbb{R}^{d}$.  Invent. Math. 220, 1 (2020), 255–393.
  * [18] Bouchut, F. Renormalized solutions to the Vlasov equation with coefficients of bounded variation.  Arch. Ration. Mech. Anal. 157, 1 (2001), 75–90.
  * [19] Bouchut, F., and Desvillettes, L. On two-dimensional Hamiltonian transport equations with continuous coefficients.  Differential Integral Equations 14, 8 (2001), 1015–1024.
  * [20] Bouchut, F., and James, F. One-dimensional transport equations with discontinuous coefficients.  Nonlinear Anal. 32, 7 (1998), 891–933.
  * [21] Bouchut, F., and James, F. Differentiability with respect to initial data for a scalar conservation law.  In Hyperbolic problems: theory, numerics, applications, Vol. I (Zürich, 1998), vol. 129 of Internat. Ser. Numer. Math. Birkhäuser, Basel, 1999, pp. 113–118.
  * [22] Bouchut, F., and James, F. Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness.  Comm. Partial Differential Equations 24, 11-12 (1999), 2173–2189.
  * [23] Bouchut, F., James, F., and Mancini, S. Uniqueness and weak stability for multi-dimensional transport equations with one-sided Lipschitz coefficient.  Ann. Sc. Norm. Super. Pisa Cl. Sci. (5) 4, 1 (2005), 1–25.
  * [24] Bressan, A. An ill posed Cauchy problem for a hyperbolic system in two space dimensions.  Rend. Sem. Mat. Univ. Padova 110 (2003), 103–117.
  * [25] Bressan, A. A lemma and a conjecture on the cost of rearrangements.  Rend. Sem. Mat. Univ. Padova 110 (2003), 97–102.
  * [26] Cardaliaguet, P., Seeger, B., and Souganidis, P. E. Mean field games with common noise and degenerate idiosyncratic noise.  arXiv:2207.10209 [math.AP].
  * [27] Cardaliaguet, P., and Souganidis, P. E. On first order mean field game systems with a common noise.  Ann. Appl. Probab. 32, 3 (2022), 2289–2326.
  * [28] Champagnat, N., and Jabin, P.-E. Well posedness in any dimension for Hamiltonian flows with non BV force terms.  Comm. Partial Differential Equations 35, 5 (2010), 786–816.
  * [29] Champagnat, N., and Jabin, P.-E. Strong solutions to stochastic differential equations with rough coefficients.  Ann. Probab. 46, 3 (2018), 1498–1541.
  * [30] Conway, E. D. Generalized solutions of linear differential equations with discontinuous coefficients and the uniqueness question for multidimensional quasilinear conservation laws.  J. Math. Anal. Appl. 18 (1967), 238–251.
  * [31] Crandall, M. G., and Ishii, H. The maximum principle for semicontinuous functions.  Differential Integral Equations 3, 6 (1990), 1001–1014.
  * [32] Crandall, M. G., Ishii, H., and Lions, P.-L. User’s guide to viscosity solutions of second order partial differential equations.  Bull. Amer. Math. Soc. (N.S.) 27, 1 (1992), 1–67.
  * [33] Crippa, G., and De Lellis, C. Estimates and regularity results for the DiPerna-Lions flow.  J. Reine Angew. Math. 616 (2008), 15–46.
  * [34] Cullen, M., and Feldman, M. Lagrangian solutions of semigeostrophic equations in physical space.  SIAM J. Math. Anal. 37, 5 (2006), 1371–1395.
  * [35] Dafermos, C. M. Generalized characteristics and the structure of solutions of hyperbolic conservation laws.  Indiana Univ. Math. J. 26, 6 (1977), 1097–1119.
  * [36] Davie, A. M. Uniqueness of solutions of stochastic differential equations.  Int. Math. Res. Not. IMRN, 24 (2007), Art. ID rnm124, 26.
  * [37] DiPerna, R. J., and Lions, P.-L. Ordinary differential equations, transport theory and Sobolev spaces.  Invent. Math. 98, 3 (1989), 511–547.
  * [38] Figalli, A. Existence and uniqueness of martingale solutions for SDEs with rough or degenerate coefficients.  J. Funct. Anal. 254, 1 (2008), 109–153.
  * [39] Filippov, A. F. Differential equations with discontinuous right-hand side.  Mat. Sb. (N.S.) 51 (93) (1960), 99–128.
  * [40] Godlewski, E., Olazabal, M., and Raviart, P.-A. On the linearization of hyperbolic systems of conservation laws. Application to stability.  In Équations aux dérivées partielles et applications. Gauthier-Villars, Éd. Sci. Méd. Elsevier, Paris, 1998, pp. 549–570.
  * [41] Godlewski, E., Olazabal, M., and Raviart, P.-A. On the linearization of systems of conservation laws for fluids at a material contact discontinuity.  J. Math. Pures Appl. (9) 78, 10 (1999), 1013–1042.
  * [42] Hauray, M. On two-dimensional Hamiltonian transport equations with $L_{\rm loc}^{p}$ coefficients.  Ann. Inst. H. Poincaré C Anal. Non Linéaire 20, 4 (2003), 625–644.
  * [43] Ishii, H. Hamilton-Jacobi equations with discontinuous Hamiltonians on arbitrary open sets.  Bull. Fac. Sci. Engrg. Chuo Univ. 28 (1985), 33–77.
  * [44] Karimghasemi, S., Müller, S., and Westdickenberg, M. Flow solutions of transport equations.  Comm. Partial Differential Equations 46, 1 (2021), 98–134.
  * [45] Keyfitz, B. L., and Kranzer, H. C. A system of nonstrictly hyperbolic conservation laws arising in elasticity theory.  Arch. Rational Mech. Anal. 72, 3 (1979/80), 219–241.
  * [46] Lasry, J.-M., and Lions, P.-L. Mean field games.  Jpn. J. Math. 2, 1 (2007), 229–260.
  * [47] Le Bris, C., and Lions, P.-L. Renormalized solutions of some transport equations with partially $W^{1,1}$ velocities and applications.  Ann. Mat. Pura Appl. (4) 183, 1 (2004), 97–130.
  * [48] Le Bris, C., and Lions, P.-L. Existence and uniqueness of solutions to Fokker-Planck type equations with irregular coefficients.  Comm. Partial Differential Equations 33, 7-9 (2008), 1272–1317.
  * [49] Le Bris, C., and Lions, P.-L. Parabolic equations with irregular data and related issues—applications to stochastic differential equations, vol. 4 of De Gruyter Series in Applied and Numerical Mathematics.  De Gruyter, Berlin, [2019] ©2019.
  * [50] LeFloch, P. An existence and uniqueness result for two nonstrictly hyperbolic systems.  In Nonlinear evolution equations that change type, vol. 27 of IMA Vol. Math. Appl. Springer, New York, 1990, pp. 126–138.
  * [51] Lions, P.-L., and Perthame, B. Remarks on Hamilton-Jacobi equations with measurable time-dependent Hamiltonians.  Nonlinear Anal. 11, 5 (1987), 613–621.
  * [52] Lions, P.-L., and Seeger, B. Linear and nonlinear transport equations with increasing velocity fields.  Forthcoming work.
  * [53] Nunziante, D. Uniqueness of viscosity solutions of fully nonlinear second order parabolic equations with discontinuous time-dependence.  Differential Integral Equations 3, 1 (1990), 77–91.
  * [54] Nunziante, D. Existence and uniqueness of unbounded viscosity solutions of parabolic equations with discontinuous time-dependence.  Nonlinear Anal. 18, 11 (1992), 1033–1062.
  * [55] Petrova, G., and Popov, B. Linear transport equations with discontinuous coefficients.  Comm. Partial Differential Equations 24, 9-10 (1999), 1849–1873.
  * [56] Petrova, G., and Popov, B. Linear transport equations with $\mu$-monotone coefficients.  J. Math. Anal. Appl. 260, 2 (2001), 307–324.
  * [57] Poupaud, F., and Rascle, M. Measure solutions to the linear multi-dimensional transport equation with non-smooth coefficients.  Comm. Partial Differential Equations 22, 1-2 (1997), 337–358.
  * [58] Suquet, P.-M. Existence et régularité des solutions des équations de la plasticité.  C. R. Acad. Sci. Paris Sér. A-B 286, 24 (1978), A1201–A1204.
  * [59] Trevisan, D. Well-posedness of multidimensional diffusion processes with weakly differentiable coefficients.  Electron. J. Probab. 21 (2016), Paper No. 22, 41.
  * [60] Veretennikov, A. J. Strong solutions and explicit formulas for solutions of stochastic integral equations.  Mat. Sb. (N.S.) 111(153), 3 (1980), 434–452, 480.
  * [61] Yamada, T., and Watanabe, S. On the uniqueness of solutions of stochastic differential equations.  J. Math. Kyoto Univ. 11 (1971), 155–167.
  * [62] Zheng, Y. X., and Majda, A. Existence of global weak solutions to one-component Vlasov-Poisson and Fokker-Planck-Poisson systems in one space dimension with measures as initial data.  Comm. Pure Appl. Math. 47, 10 (1994), 1365–1401.