Option prices under liquidity risk as weak solutions of semilinear diffusion equations

Prices of financial options in a market with liquidity risk are shown to be weak solutions of a class of semilinear parabolic partial differential equations with nonnegative characteristic form. We prove the existence and uniqueness of such solutions, and then show the solutions correspond to option prices as defined in terms of replication in a probabilistic setup. We obtain an asymptotic representation of the price and the hedging strategy as a liquidity parameter converges to zero.


Introduction
Nowadays, risk cannot be efficiently managed without taking into account liquidity risk. An important aspect of risk management for financial institutions is to understand the effect of liquidity on the pricing and hedging of derivative securities. The definition of liquidity depends on the market structure and the financial or economic questions being studied. For instance, in a limit order book market, one can measure liquidity by quantifying the number of shares being offered at each price, or the resilience of the order book after a large trade. On the other hand, to study liquidity in a dealers market, one must typically consider some notion of information flow, noise trading and utility functions for the market makers.
In this paper, we consider liquidity in terms of the depth of the market, namely the impact that trades have on prices. We give a framework for the problem of option pricing in a large trader model with a partial differential equation (PDE) perspective. We analyse the limit as the market becomes infinitely liquid (the impact of a trade becomes infinitely small) and provide first-order asymptotics of the option price and the hedging strategy in a liquidity parameter. For the economic model we follow Jarrow and Roch [29] and consider a variation based on Brownian motion of the larger trader model of Bank and Baum [3] in which the dynamics of the price process depends specifically on the current holdings of the investor. The price process depends on economic variables which follow a degenerate diffusion process. Unlike [3] however we do not specifically assume the existence of a local martingale measure for all primal processes (Assumption 3 in [3]). The consequence is a nonlinear term in the wealth equation and a liquidity premium in option prices. Our analysis also has some common features with the liquidity cost model of [10], and in particular the Taylor expansion of the super-hedging cost of [39] who use the PDE characterisation of [11]. The main difference however is that in these models the price impact is momentary, whereas we regard it as a longer-lasting phenomenon.
In mathematical terms we derive a semilinear parabolic PDE on a bounded domain in R n given as where the second-order partial differential operator is in divergence form for a matrix σ and a vector b. We assume that the quadratic form defined by σσ , i.e. the characteristic form, is nonnegative. Moreover, the error term F , representing the liquidity cost, is quadratic in σ Du. The characteristic form being nonnegative as opposed to bounded below by a positive number (i.e., uniform ellipticity of the operator L) is a type of degeneracy. Using variational methods we show existence and uniqueness of weak solutions of (1). We work with the concept of weak solutions for the simple reason that it allows us to obtain information about the growth of the gradient σ Du in an L 2 -space. Many authors have approached the problem of existence and uniqueness of PDEs, and the relation to backward stochastic differential equations (BSDEs) with the theory of viscosity solutions [13,23], but this type of solution is not differentiable in general.
A clear advantage of the weak solution approach is that we obtain the (economically crucial) hedging strategy as the gradient of the solution, which allows us to show that the hedging strategy converges to the hedging strategy in a frictionless setting as the market becomes more liquid. Finally, we investigate the regularity of the option price u and σ Du and show that these functions are Hölder continuous on certain subdomains.
We briefly summarise the related PDE literature distinguishing between analytic and stochastic perspectives. Due to the vastness of the literature we mainly concentrate on semilinear PDEs with quadratic growth in the gradient noting that other types of growth have also received considerable attention.
The reader is invited to consult the cited papers for further references. We remark that the study of the well-posedness of uniformly parabolic semilinear problems is standard, cf. [2,33,35], and for the case of nonnegative characteristic form we refer to [36]. To paraphrase [4], this is, however, not the kind of degeneracy that has received wide attention compared to the type found in quasilinear equations where the coefficients of the second order operator depend on u and Du.
To the best of our knowledge there is no treatment of parabolic PDEs of nonnegative characteristic form with error terms of quadratic growth. In [28], the author treats quasilinear PDEs of a special structure and growth conditions which do not cover (1). Closely related to our problem is [44] which which treats existence, uniqueness and regularity of weak solutions of PDEs of type (1), however under a uniform ellipticity assumption. More recent attention focused on the quasilinear case under a coerciveness assumption, cf. [7,8,21]. A quasilinear initial-boundary-value problem also in an economic context using the assumption of uniformly coercive operators was studied using weak solutions in [16].
From a stochastic perspective, the connection between BSDEs and semilinear partial differential equations goes back to [38] and was further exploited in terms of classical and viscosity solutions in [37]. The case of quadratic growth in the gradient was considered in [30] in two ways: the elliptic Dirichlet problem for a uniformly elliptic operator on a bounded domain was treated using weak solutions and the Cauchy problem on Euclidean space using viscosity solutions. Important for our paper is the seminal [5] where the link between BSDEs and weak solutions of the associated PDEs is explored. Here, the generator is assumed to be at most linear in the gradient Du and the authors consider the Cauchy problem on Euclidean space.
In the more recent literature, a generator F (τ, x, u, Du) with at most quadratic growth in u was treated in [40] using weak solutions on Euclidean space. Linear degenerate initial-boundary value problems on unbounded domains (also motivated by financial mathematics) are considered in [19,20] and further papers by the same authors.
This paper is organised as follows. The next section introduces the notation. Section 3 develops the probabilistic formulation for our liquidity models.
The key results of the PDE analysis are stated and interpreted in Sect. 4 with detailed proofs contained in Sect. 5.

Notation
We introduce certain spaces and pieces of notation. The reader is referred to [1,18] for further details and proofs.
Let Ω ⊂ R n be a bounded domain, denote by vol(Ω) the Lebesgue measure of Ω and define Q T = [0, T ] × Ω to be the parabolic cylinder. We denote by C(Ω) the linear space of bounded continuous functions on Ω and by C c (Ω) the linear space of continuous functions that have compact support. Both are metric spaces under the usual sup norm ||f || ∞ = sup x∈Ω |f (x)|, analogously for C(Q T ) under the norm ||f || ∞ = sup (t,x)∈QT |f (t, x)|.
For k ∈ N we let C k (Ω) be the space of k times continuously differentiable functions with norm ||f 1 , . . . , D αn n ) for α = (α 1 , . . . , α n ) ∈ N n 0 a multi-index, |α| = α 1 + · · · + α n , and D i = ∂ i = ∂/∂ i denotes differentiation with respect to the ith coordinate. By D = (D 1 1 , . . . , D 1 n ) we denote the gradient operator. The analogous notation applies for functions of time and space where we have the function spaces Let δ be a metric on Ω. The Hölder continuous functions with exponent α ∈ (0, 1) are given by Typically one uses the metric δ(x, y) = |x − y| with the Euclidean distance between x and y. The analogous definition holds on Q T where we mention the special metric the so-called parabolic metric.
The space C ∞ c (Ω) is the space of test functions, i.e. smooth functions of compact support. The Lebesgue spaces L p (Ω) are defined as the set of equivalence classes of measurable real-valued functions f on Ω such that ||f || p = Ω |f (x)| p dx is finite; similarly for L p (Q T ). By L ∞ (Ω) we denote the space of measurable functions with the essential supremum norm.
Given f, g ∈ L 1 (Ω) and α a multi-index, we say g = D α f in a weak sense (αth weak partial derivative) if Ω fD α ϕdx = (−1) |α| Ω gϕdx for all ϕ ∈ C ∞ c (Ω). The Sobolev spaces W k,p (Ω) for k ≥ 1 and 1 ≤ p < ∞ are then defined as the set of measurable functions f with weak derivatives D α f up to order k such that the norm ||f || p . For X a Banach space of real-valued functions on Ω we define the space C([0, T ]; X) to be the set of functions f : We denote the inner product in L 2 (Ω) by (·, ·). The inner product in L 2 (Q T ) (and by abuse of notation also in Given a bounded domain Ω with C 1 -boundary ∂Ω we have the integration by parts formula with ν = (ν 1 , . . . , ν n ) the outward pointing unit normal field and dS a surface element.

NoDEA
Option prices under liquidity risk as weak solutions Page 5 of 32 12 Given a second-order linear differential operator with suitable coefficient functions a ij , b i , we say that L(τ ) is uniformly elliptic on Ω if the quadratic form given on for all τ ∈ [0, T ] and (x, ξ) ∈ Ω × R n , where | · | is the Euclidean norm on R n .

Probabilistic setup
We consider an investor who incurs liquidity costs due to the trading impact on prices of a risky asset. We present the price impacts model and define the replication problem for which the solution of the PDE will be shown to be associated to option prices.

The price impact model
Following [29], we consider a multi-asset variation of the large trader model of [3], which we base on Brownian motion. The interest rate is constant, and for simplicity we only consider discounted prices. We consider a market that consists of d traded risky assets. Economic variables. We first define an n-dimensional process (X t ) t≥0 that represents all economically relevant variables, i.e. fundamentals prices, volatilities, interest rates, market liquidity, etc. To this end let W be a ddimensional Brownian motion, and F = (F t ) t≥0 its filtration. The process (X t ) is the unique strong solution of with β and σ Lipschitz continuous in the second argument uniformly t. Here, β is R n -valued and σ takes values in the n × d-matrices over R.
Price of risky assets. We assume there is a family of d-dimensional stochastic processes S(θ) = (S(t, X t ; θ)) t≥0,θ∈R d for which each component of S(t, X t ; θ) is interpreted as the price of a traded risky asset at time t when the investor holds a constant position of θ ∈ R d units of these assets. The i-th component of θ, denoted θ i , gives the position in the i-th asset. In this sense, holding one of the d assets may have a price impact on any of the d assets. A negative value for θ i represents a short sale in this asset. We often write S t (θ) for S(t, X t ; θ) to simplify the notation.
To derive the dynamics of the prices, we suppose that S(·; θ) : [0, T ] × R n → R d is a deterministic function that is continuously differentiable with respect to t and twice continuously differentiable with respect to x. By Itô's Formula, the family of processes can be expressed with the following stochastic differential equations We assume there is a local martingale measure Q for the unaffected price processes (the price processes observed if θ ≡ 0): This assumption rules out arbitrage opportunities for small traders (traders who do not have an impact on prices), when the large trader does not trade. In terms of B, X takes the representation , we do not assume that all S(θ) are local martingales under Q. Indeed, for all θ, we assume there is a local martingale measure Q(θ) for the price processes S i (θ): is the risk premium associated to the i-th Brownian motion risk W i when the investor holds a position θ. The above equation implies that the large trader has a direct impact on the risk premia of the traded assets.
The existence of η(·; θ) is justified by the fact that when θ is kept constant by the large trader, small traders obtain the price S i t (θ), so the existence of the equivalent local martingale measure Q(θ) rules out arbitrage opportunities for small traders, cf. Hypothesis 2 (NFLVR Infinitesimal Traders) in [29]. With these measure changes, S i (θ) is represented as

Example 1.
A simple example is the Bachelier model for θ = 0 under Q(0): with σ a d × d matrix, and Option prices under liquidity risk as weak solutions Page 7 of 32 12 For this model,σ ji (t, X t ; θ)S i t (θ) = σ ji , and ψ t (θ) = 2λσθ. A well-known empirical feature of asset prices is that the risk premium depends on the volatility (cf. [24]). In this simple model, the large trader's impact on the risk premium is proportional to the volatility. The i-th component of the vector σθ gives the large trader's exposure to the i-th Brownian risk W i so that the change of drift associated to the i-th component of W is proportional to this exposure.
A first trade at time t makes the price process S i jump by 2λ(σ i ) σθ(T − t). Also, at time T , all processes S i (θ) converge back to S i (0). Liquidity costs. We follow [3] who define the asymptotic liquidation proceeds from a position θ at time t for the single asset setup by considering that the asset is liquidated in infinitesimal packets, infinitely fast. In other words, the liquidation value of an asset at time t is given by the integral In [29], this definition is extended to d-dimensional trading strategies. The definition of the asymptotic liquidation proceeds L t (θ) then involves a curvilinear integral from θ to the d-dimensional vector 0. To simplify the treatment, we adopt the convention that each asset is purchased (resp. liquidated) one by one, in the (resp. reverse) order given by their index i. Shares of the assets are liquidated in infinitesimally small packets, starting from asset d, down to asset 1. However, since prices are a function of the current holdings, the liquidation of the i-th asset is executed while still holding the first i − 1 assets. Hence, the price obtained during this liquidation is S i t (Θ i (y)) and the liquidation value of asset i is with Θ i (y) = (θ 1 , θ 2 , . . . , θ i−1 , y, 0, . . . , 0), (y ≤ θ i ). As such, for a position θ = (θ i ) i≤d to liquidate, we define the asymptotic liquidation proceeds as Following the same logic, the cost of building a position θ (a negative quantity) is defined as by considering that assets are purchased in the order of their index. A simple change of variable shows that this is simply equal to −L t (θ). More generally, the proceeds obtained from changing position from θ to ϑ is given by L t (θ) − L t (ϑ) at time t.
By Theorem 2.2 of [45] (Fubini's Theorem for stochastic integrals) we can write the asymptotic liquidation proceeds as with and the components of the vector-valued process Σ are given by

Trading strategies, and wealth processes
We now define the wealth processes associated to a self-financing trading strategy. As opposed to [3] who go through considerable details to define the wealth process, the representation for L in (5) allows us to take an easier route. This construction is taken from [29] and works for any value of d ≥ 1.
We start by defining the notion of self-financing for simple trading strategies. Let θ t = i≥1 ξ i 1 {τi≤t} be a simple strategy in which (τ i ) i≥1 is an increasing sequence of stopping times and ξ i ∈ F τi . The vector θ t denotes the number of shares owned by the investor in each risky asset at time t. We define the wealth process Π directly in terms of L. As shown in [3], this allows the investor to minimize transaction costs associated to the liquidity of assets. More precisely, at time t, the wealth is given by Option prices under liquidity risk as weak solutions Page 9 of 32 12 By the fact that θ τi−1 = θ τi− for i > 1 and θ 0 = 0, the sum becomes using (5) and the fact that θ s = θ τi−1 for τ i−1 ≤ s < τ i . Because these two integrals are well-defined for more general processes θ, we can extend the definition of wealth processes to more general trading strategies: for T > 0. The wealth process associated to θ is then given by

Replication of contingent claims
The main problem in which we are interested is the pricing of contingent claims in the context of liquidity risk and price impacts. We proceed by calculating the replication cost of contingent claims payoffs. Let T > 0 denote the maturity of an option. Ifĥ denotes its payoff function at time T , then the replication problem consists in finding a trading strategy θ that sets the portfolio wealth at time T equal toĥ(S T ): Remark 5. If an investor owns the option and wants to hedge away the risk, he needs to solve: which corresponds to (7) withĥ replaced with −ĥ.
Of course, the expressionĥ(S T ) is ambiguous here as it depends on θ T . In order to avoid price manipulations, however, it is often assumed that the investor liquidates his position at time T , so that the observed price at maturity is S(T, X T ; 0) and the associated payoff isĥ(S(T, X T ; 0)), or that θ T is replaced by an approximation Δ(T, X T ). In general, the option payoff can thus be represented as h(X T ), for some deterministic function h.
We can make this a dynamic problem by making the following definition: in which θ is called the replication (or hedging) strategy.
Equation (8) is a type of BSDE for which the pair (Π, θ) is a solution.
We make the following standing assumption on Σ: Under this assumption, for each t ≤ T, z ∈ R d , x ∈ R n we can find θ 0 ∈ R d such that Σ(t, x; θ 0 ) = z. Accordingly, we can find a measurable function Φ (see Lemma 1 of [29]) such that and Φ will play a key role in the PDE analysis.
In Example 1, Assumption 7 is satisfied if and only if σ is invertible. In this case, Φ(z) = λ|z| 2 .
The quadratic form for Φ can be obtained for much more general models than Example 1. In fact, for any choice ofσ and μ it suffices to take in whichf kj =f jk for all 1 ≤ k, j ≤ d. In this case, consider the change of variable z j = Σ j (t, x; θ) so that is a conservative vector field, it follows that from the fact that ∇Φ = φ. From an economic point of view, the representation of ψ as in Eq. 10 means that the required market premium for Brownian motion W j is a linear combination of the volatility structure Σ(t, x; ·), thus giving a multivariate generalisation of the well-known empirical feature of asset prices that risk premia depend on the volatility (cf. [24]). We make this a standing assumption.

Hypothesis 8.
We assume that In practically relevant situations, one could directly specify the parameters in the dynamics of X and the function Φ, (without the need to specifyσ and μ) reflecting the properties of a given financial market (economic variables and liquidity costs associated to trading).
Hypothesis 7 implies that the existence of a solution (Π, Z) of the BSDE (8) is equivalent to the existence of a solution (Y, Z) of Since this is not a linear equation in Z, the replication cost of two units of h is not twice the replication cost of one unit. In order to emphasise the dependence on this nonlinear term and to study the asymptotic representation of Y and Z when liquidity costs are small, we introduce the parameter λ > 0 in the previous equation: By Theorem 2 of [9], the BSDE (11) has a solution when z → Φ(t, x; z) is continuous for all t, x and has at most quadratic growth in z, i.e. when there are constants C 0 , C 1 ≥ 0 such that |Φ(t, x; z)| ≤ C 0 + C 1 |z| 2 uniformly in t, x and E Q e C1h(XT ) < ∞.
To this general replication problem we naturally associate a nonlinear PDE. We prove the existence and uniqueness of the solution of this PDE in a Sobolev space, and show in a second step that it also gives the solution of the BSDE. The clear advantage of this approach is that we explicitly obtain the replication strategy as the derivative of the option price with respect to the underlying. This also allows us to show that the replication strategy and the option price converges to the replication strategy and the option price in a frictionless setting when liquidity costs are small (λ is small).
In the setting of Example 1, Φ is given by Φ(z) = λ|z| 2 . Consider the variable χ = e 2λY . By Itô's Formula and Eq. (11), In other words, χ is a martingale, and χ t = E Q exp(2λh(X T )) F t . Consequently, Y can be represented as The solution Y can therefore be represented explicitly in terms of h(X T ) in this special case.

Bounded domains and the PDE formulation
Our goal is to obtain the solution of the above replication problem in terms of an associated PDE, and study the analytical properties of the solution. The PDE results below are valid for a bounded domain Ω with C 2 -boundary ∂Ω. This is a common assumption in the PDE literature, and our results below cannot be extended to the case of domains of infinite volume. In practice, the solution of the PDE gives a good approximation of the problem as one can take Ω as large as needed. Indeed, suppose that the initial condition of the process X at time t is given by X t = x a.s. Introducing the stopping time τ If Ω m is a domain that contains the ball of radius m centered at 0 and sup t≤T Υ(t, The solution of (12) thus satisfies the approximation replication definition of [10]. In this sense, one can make the mean square hedging error arbitrarily close to zero, and that should be quite satisfactory for any practical implementation, in particular when using finite difference methods. A simple specification for the function H is of course obtained by setting g = 0. Note that [3] defined the related notion of approximately attainable contingent claim with a similar line of reasoning.
After time reversion τ = T − t, the PDE associated to (12) is with Here, Ω is the domain of X, and b = (b 1 , . . . , b n ) and σ = (σ ij ) 1≤i≤n,1≤j≤d are the drift vector and diffusion matrix of the process X. When λ = 0, the market is perfectly liquid (the price impact of a trade is zero). In this case, option prices (solutions of (12) or (13)) are Q-martingales and are related by which is the classical Feynman-Kac result that tells us that the solution of (13) can be written as a conditional expectation of the terminal condition.

NoDEA
Option prices under liquidity risk as weak solutions Page 13 of 32 12

Analysis of the PDE
We address existence and uniqueness of solutions to the PDE (13), first-order asymptotics in λ and the Hölder continuity of solution (and their gradients) on certain subdomains of Ω. We consider the slightly more general initial-boundary value problem where L(τ ) is written in divergence form with detailed assumptions to be given below. We understand u(τ, x) as a Banach-space-valued function of τ so that we mostly drop any x-dependence. The expression σ (τ, x)Du(τ, x) for a function u will always be treated as a single entity which we abbreviate as σ Du. This is the σ -gradient in the sense of Chapter 4.1 of [28].

Hypothesis 9.
We make the following assumptions.
(i) Domain: the set Ω is open, bounded and has a boundary of class C 2 .
(ii) Coefficients: The n × d-matrix σ has components σ ij which belong to C 0,2 (Q T ) and the components of the drift vector b = (b 1 , . . . , b n ) belong to C 0,1 (Q T ). The quadratic form defined by the square matrix σσ is not assumed to be positive definite. (iii) Generator: The function F takes the form Since the characteristic form of L may be zero, we seek solutions in a Sobolev space based on the σ -gradient. This space is defined as Standard arguments (cf., Theorem 3.3 of [1]) show that W 1,2 σ (Ω) is a Banach space. Also let W 1,2 σ,0 (Ω) be the closure of C ∞ c (Ω) in W 1,2 σ (Ω). The solution of the PDE (15) will then be in the Hilbert space L 2 ([0, T ]; W 1,2 σ (Ω)) and the boundary condition will be interpreted in terms of L 2 ([0, T ]; W 1,2 σ,0 (Ω)).
Also recall the notion of a weak (or generalised) solution.

NoDEA
Option prices under liquidity risk as weak solutions Page 15 of 32 12 Here B is a function of τ given by the bilinear form The initial-boundary condition is interpreted as u(0) = g in L 2 (Ω) and u−H ∈ L 2 ([0, T ]; W 1,2 σ,0 (Ω)). Remark 13. The interpretation of the boundary condition looks unusual since one typically phrases this in terms of the trace operator. However, the space W 1,2 σ (Ω) does not possess a straightforward trace mapping into L 2 (∂Ω). For the classical trace W 1,2 (Ω) → L 2 (∂Ω), the kernel of this map is given precisely by W 1,2 0 (Ω) and it is this analogy that we exploit. Note that even if the weak solution u were continuous, we cannot generally expect u = g in a pointwise sense on the boundary. This is for two reasons. First, we do not assume that σ is a square matrix. Second, σ could be zero on the boundary so that near this part of ∂Ω we cannot control the gradient Du.
In the case of square σ, the existence of suitable trace mappings on the energy spaces with degenerate weighting is discussed in detail in Chapter 4.2 of [28]. We also refer to [22] and to Chapter I.1 of [36] for a detailed discussion on how to decompose the boundary into singular and regular parts based on the Fichera function.

A Feynman-Kac-type theorem
We first show the usefulness of studying this PDE: weak solutions of (15) can be used to obtain solutions of the BSDE (12) so that studying the PDE solutions yields valuable information about option prices and their gradients. This is expressed in a Feynman-Kac-type theorem the proof of which is nontrivial as we do not deal with classical PDE solutions, i.e. twice continuously differentiable functions. For ease of notation we suppress the initial condition X t = x a.s. from the statement of the theorem. Theorem 14. Assume Hypothesis 9 and let u ∈ L 2 ([0, T ]; W 1,2 σ (Ω)) be a weak solution of the initial-boundary-value problem (15).
The proof of the theorem is strongly intertwined with the existence and uniqueness result for the PDE so that it is contained in Sect. 5.5.

Existence and uniqueness
We now state the existence and uniqueness of solutions to the aforementioned class of semilinear PDEs on bounded domains. In economic terms the result shows that there is a unique option price in an L 2 -space whose σ -gradient (on which the hedging of the contingent claim is based) also lives in an L 2 -space. The latter assertion cannot be obtained through viscosity solutions. The proof relies on the Schauder fixed point theorem. In numerical implementations, a Galerkin scheme could be used to solve a regularised PDE and pass to the weak limit.

Liquidity asymptotics of the solution
In this section we describe the continuity of the PDE solutions with respect to the parameter λ leading to natural asymptotics. From an economic point of view, we study the marginal properties of prices by considering the limit as the assets become more liquid.
Theorem 16. Assume Hypothesis 9 and let u (λ) be the unique weak solution of the PDE (15) for a suitably small λ ≥ 0.

(i) Then as
The first assertion expresses the continuity from above at 0 of the derivative prices as a function of the liquidity parameter: as the market becomes more liquid, the derivative price continuously approaches the Black-Scholes price (14) of the derivative when price impact is absent.
The second assertion also makes precise a formal perturbation approach in powers of λ. The function v (λ) gives the additional liquidity cost in the option price u (λ) per unit of λ, i.e. the marginal liquidity cost of the option.

Regularity of the weak solution
In practical implementations one is naturally interested in the regularity of the weak solution of the PDE (15). The more regular the PDE solution, the better is the convergence of such a numerical scheme.
Hölder regularity of parabolic equations and systems is a topic wellcovered in the literature. For quasilinear equations with at most quadratic growth in Du we refer to Chapter V.1 of [33] whose Theorem 1.1 establishes Hölder continuity given a smallness condition, cf. also [12] which shows Hölder continuity once Du is in an L q -space. Other types of growth in Du are covered extensively in [15]. The corresponding results for semilinear parabolic system with quadratic growth in the gradient can be found in [26,27,43,44] to name just a few references. A broader overview is given in [31].
where ||A|| op denotes the operator norm of a real d × d-matrix A acting on R d . Then the weak solution u and its gradient Du are Hölder continuous on Ω with respect to the parabolic metric δ defined in (2) for some Hölder exponent α ∈ (0, 1) which depends on the data of the PDE and on ν.

Proofs of the key results
This section is more technical as it contains the proofs of the key results. Each subsection contains one building block of the proofs, blocks 2.-5. depend on the estimates derived in 1. For ease of presentation we only consider time-independent generators of the diagonal form F : Ω × R d × R d → R is given as

Existence and uniqueness
where we set γ = max 1≤i≤d ||f i || ∞ . This is no restriction and the general case follows similarly in each instance.

Existence and uniqueness result with zero boundary conditions
This section contains the proof of Theorem 15 in the special case when the solution is required to vanish on the boundary. The strategy of the proof is the same as in [5] with two nontrivial complications: the quadratic nonlinearity and the boundary.
To solve the PDE we interpret it in a variational sense and use energy methods based on L ∞ -a-priori estimates obtained via a classical maximum principle. Since the differential operator L is independent of the solution u, we do not need more abstract methods to treat quasilinear and fully nonlinear equations, cf. [41,47]. So our exposition is self-contained and easily accessible.
The diffusion degeneracy requires a regularisation of the equation by the method of vanishing viscosity (elliptic regularisation), cf. for example [6,34,36,42] and the weak convergence of the corresponding solutions [17,18]. This approach is also suggested by the connection with stochastic analysis. We want to establish the link between a class of evolution equations and BSDEs. This is achieved by suitably approximating the PDE problem so that the classical Feynman-Kac theorem can be used. Proof. The idea (cf. Chapter 14 of [5]) is first to obtain a priori estimates of u in L ∞ (Q T ) for the simplified problem for some fixed v ∈ L 2 ([0, T ]; W 1,2 σ (Ω)). These estimates then allow the application of weak convergence methods to construct a solution of (18). Denoting the solution of (18) by u = A [v] to highlight the dependence on v we define a nonlinear operator A. This turns out to be a compact operator preserving a suitable subset of a Hilbert space so that the fixed point of A guaranteed by the Schauder fixed point theorem is the desired solution to the original PDE. The precise argument is as follows: (i) Let v ∈ L 2 ([0, T ]; W 1,2 σ (Ω)) and approximate σ Dv in the L 2 -norm by a sequence of smooth functionsṽ (Lemma 19) indexed by > 0. Consider the corresponding uniformly elliptic PDE Here, Δ stands for the Laplace operator on R n . By standard results (e.g., Theorem 4 in Chapter 9 of [25]), the PDE (19) has a unique classical solution u . First of all note that using mollifiers, any non-continuous h can be approximated arbitrarily closely in the L 2 (Ω)-norm by a sequence of smooth functions h . By construction, these function satisfy ||h || ∞ ≤ ||h|| ∞ . This allows us to extend all proofs to the case of h ∈ L ∞ (Ω). The same argument applies to f . Expressing u using a Duhamel representation based on a fundamental solution shows that we can allow for h and f measurable. (ii) Note that due to the linearity of F in z we can "pull" the term in (iii) Using these a priori estimates and invoking a compactness argument we can extract a weakly converging sequence u ∈ L 2 ([0, T ]; W 1,2 σ (Ω)) ∩ L ∞ (Q T ) by Lemma 22 and Corollary 24. This sequence also converges strongly to a limit u ∈ L 2 ([0, T ]; W 1,2 σ (Ω)) as shown in Lemma 23. The key point here is that ||divb|| L ∞ (QT ) < ∞ which restricts the choice coefficients b i and σ ij . Moreover, this function is regular in the time variable in the sense that u ∈ C([0, T ]; L 2 (Ω)). (iv) So far, the solution u depends on v. This allows us to define a continuous nonlinear map A : v → u defined as the strong limit of the u wheñ v → σ Dv. We show that that A is well-defined (Lemma 25) and for λ < 1/2γm preserves the space Z (Lemma 26), where with R sufficiently large depending on the data of the PDE. Since A is continuous and compact (Lemma 25), it has a fixed point in Z by the Schauder fixed point theorem, cf. for example Theorem 2.A of [46]. This fixed point is a weak solution of the PDE (Lemma 27).
(v) The PDE has at most one solution: the arguments of Appendix A.2 in [5] go through in our situation. (vi) It remains to show that for f = 0 the existence claim holds for any λ ≥ 0. Suppose that λ < 1/2γ||h|| ∞ and choose μ > 1. Since λ < 1/2γ||h|| ∞ < 1/2γ|| 1 μ h|| ∞ , by the above arguments the PDE also has a unique weak solution u. The homogeneity of F implies that the function v = μu is a weak solution of Since μ > 1 was arbitrary, the problem (15) has a solution for any λ.
This completes the proof.
For an element in L 2 ([0, T ]; W 1,2 σ (Ω)) we need to approximate the weak derivative σ Dv in a controlled way. We shall also employ the following classical maximum principle.  As a final preliminary result we note a useful bound on the inner product of F with L ∞ -functions. This follows from Cauchy's inequality.
Using Gronwall's lemma we obtain estimates for λ sufficiently small.
The existence of the weak limits u and σ Du is related to the weak closure of the σ -gradient operator, cf. the discussion preceding Lemma 1.1 of Chapter 4.1 in [28].
1. Let , > 0 and suppose thatṽ ,ṽ ∈ L 2 (Q T ) ⊗ R d as in Lemma 19. As before, let u , u be the corresponding solutions of (19). We find The tricky term in this expression is .
. (22) 2. Via Gronwall's inequality we see that u converges strongly to u. This is due to two reasons. First, after integrating (22) with respect to τ , the term 2 Δu , u + 2 Δu , u tends to zero as , → 0 since 2 Δu converges weakly to zero in L 2 (Q T ) (the arguments of the proof of Lemma 14.8 of [5] go through, cf. also the proof of Theorem 5 of Chapter 7.1.3 in [18]). Second, the inner product in F tends to zero using the bilinearity of F , the L 2 (Q T ) ⊗ R dconvergence ofṽ and the weak convergence of σ Du . Moreover, the bound on ||(u − u )(τ )|| from Gronwall's inequality can be chosen independently of τ so we have convergence in C([0, T ]; L 2 (Ω)).
3. As regards the convergence of σ Du − σ Du in L 2 (Q T ) ⊗ R d , we see from (22) that ||σ Du − σ Du || 2 L 2 (QT )⊗R d also tends to 0 as , → 0. The following corollary summarises important bounds on u. Corollary 24. In the above notation the following holds. Proof. (i) As the sequence u converges to u in L 2 (Q T ), there is a subsequence that converges a.e. to u. Since the subsequence is uniformly bounded in L ∞ (Q T ) by ||h|| ∞ + ||f || ∞ , it follows that u ∈ L ∞ (Q T ). The function u depends on the given v and we express this correspondence by defining a nonlinear map A setting u = A [v]. We have seen that this map A acts on L 2 ([0, T ]; W 1,2 σ (Ω)) and we collect key properties of A. (ii) Let v k be a convergent sequence of functions in L 2 ([0, T ]; W 1,2 σ (Ω)) with limit v. Define u k = A[v k ] and u = A [v]. We know from Corollary 24 that the sequence (u k ) is bounded in norm uniformly in k. So there is a weakly convergent subsequence (u kj ) that converges to some w. By the definition of weak solutions in Definition 12 we see that w = A [v]. Using approximating sequences u kj → u kj and u → u the argument in the proof of Lemma 23 shows that u kj → u strongly. A contradiction argument shows that we must have u k → u for the whole sequence.
(iii) The compactness claim follows by a similar argument also exploiting the uniform boundedness of A[v k ] and the consequent existence of a strongly convergent subsequence.
In order to apply the Schauder fixed point theorem we must first show that A preserves a smaller set within L 2 ([0, T ]; W 1,2 σ (Ω)).

Proof of the existence and uniqueness for nonzero boundary conditions
We extend the existence and uniqueness result to nonzero boundary conditions exploting the compatibility expressed by the function H. Proposition 28. The assertion of Proposition 18 also holds when the boundary condition g is nonzero and satisfies Hypothesis 9.
Proof. The idea of the proof is to reduce the more general boundary value problem to a problem with zero boundary conditions, cf. Remark 6.4 of [6]. This exploits the fact that F is given by a quadratic form in Du.
Suppose that u is a weak solution of (15)  We must show that this PDE can be covered by our existence and uniqueness result of Proposition 18. We assumed in Hypothesis 9 that ∂ τ H and ∂ i ∂ j H all belong to C(Q T ). The operator L can be brought into the form where the coefficients b i are such that ∂ j b i ∈ C 0,1 ([0, T ]×Ω) by the assumptions on H. The functions b i are obtained from the b i and the second summand in (27).
The converse also holds, i.e. when v solves (26), then u = v + H solves (15) with the correct boundary conditions. We have chosen all assumptions so that (26) has a unique weak solution.