Christophe Raynaud

Eliott Norbert-Agaisse

Brian R James

Odile Eisenstein

31 P Chemical Shifts in Ru(II) Phosphine Complexes. A Computational Study of the Influence of the Coordination Sphere

HAL is

INTRODUCTION

One of our groups first reported in 1993 on the square pyramidal complexes trans-RuCl2(PN)(PR3) (R = Ph, p-tolyl; PN = o-diphenylphosphino-N,N'-dimethylaniline), and their ability to bind H2 and N2 to form, as final products, the six-coordinate cis-RuCl2(PN)(PR3)L (L = H2, N2) complexes. [START_REF] Mudalige | Molecular-Hydrogen, -Nitrogen and Monohydride Derivatives of the Structurally Characterized Dichloro(odiphenylphosphino-N,N-dimethylaniline) [tris(p-tolyl)phosphine]ruthenium(II) Complex[END_REF] Subsequently, the other relatively small molecules H2S, 2,3 N2O, [START_REF] Pamplin | The Nitrous Oxide complex RuCl2(h 1 -N2O)(PN)(PPh3) (PN= [o-(N,N-Dimethylamino)phenyl] diphenylphosphine); Low Temperature Conversion of N2O to N2 and O2[END_REF] thiols, [START_REF] Ma | Ruthenium(II) Thiol and H2S Complexes: Synthesis, Characterization and Thermodynamic Properties[END_REF] CO, [START_REF] Ma | Binding of CO and NH3 at a five-coordinate Ru(II) centre in the solid state and in solution[END_REF] NH3, [START_REF] Ma | Binding of CO and NH3 at a five-coordinate Ru(II) centre in the solid state and in solution[END_REF] and HCCPh, which is converted to coordinated C=CHPh, [START_REF] Ma | Reactions of phenylacetylene and p-tolylacetylene with a five-coordinate Ru II complex[END_REF] reacted similarly to form related six-coordinate complexes, in which the chloride ligands are either cis or trans. Scheme 1 shows the cis and trans products for these systems, although the L = H2O and alcohol species were made via reactions with 1 RuCl2(PN)(PR3) that was formed in situ from RuCl3(PR3)3 complexes; [START_REF] Ma | Reversible binding of water, methanol, and ethanol to a five-coordinate ruthenium (II) complex[END_REF] the most recent report describes a corresponding cationic cis-NO + species that was formed by reaction of fac-RuCl3(NO) (PN) with PPh3. [START_REF] Da Silva | Studies on the Ru(II) monocationic complexes [RuCl2)(NO)(PN)(PPh3)]PF6, where PN = odiphenylphosphino-N,N-dimethylaniline, and R = Ph and p-X-C6H4 (X = OMe, Me, F)[END_REF] The two Scheme 1. Schematic reactivity of trans-RuCl2(PAN)(PR3) toward small molecules (L) to form isomeric six-coordinate complexes. Note: the cis-C=CHPh complex is formed by addition of HCCPh to the five-coordinate complex, and the NO + complex by reaction of PPh3 with fac-RuCl3(NO)(PN) P-atoms in the complexes are labelled PA and P (Scheme 1), as in the associated [START_REF] Jameson | Absolute Shielding Scale for 31 P from Gas phase NMR Studies[END_REF] P{ 1 H}-NMR spectra, which show the expected doublet of doublets pattern. Crystals suitable for X-ray studies have been grown for many complexes, and an empirical inverse correlation between the [START_REF] Jameson | Absolute Shielding Scale for 31 P from Gas phase NMR Studies[END_REF] PA chemical shift and the RuPA bond length has been established over the respective ranges of about 30 80 ppm and 2.392.17 Å (Figure 1). [START_REF] Da Silva | Studies on the Ru(II) monocationic complexes [RuCl2)(NO)(PN)(PPh3)]PF6, where PN = odiphenylphosphino-N,N-dimethylaniline, and R = Ph and p-X-C6H4 (X = OMe, Me, F)[END_REF] Of note, the complexes with different X groups (Cl or Br) and phosphines (PPh3, P(p-tolyl)3), or P(p-FC6H4) are present on this correlation line. Analogous trends exist for a more diverse range of Ru II -complexes containing PPh3 [START_REF] Jessop | Hydrido Thiolato and Thiolato Complexes of Ruthenium(II) Carbonyl phosphines[END_REF] and Ph2P(CH2)4PPh2 ligands, [START_REF] Queiroz | The reactivity of five-coordinate Ru(II) (1,4bis(diphenylphosphino)butane) complexes with the N-donor ligands: ammonia, pyridine, 4substituted pyridines, 2,2′-bipyridine, bis(o-pyridyl)amine, 1,10-phenanthroline, 4,7diphenylphenanthroline and ethylenediamine[END_REF][START_REF] Macfarlane | Characterization of Five-Coordinate Ruthenium(II) Phosphine Complexes by X-ray Diffraction and Solid-State 31 P CP/ MAS NMR Studies and Their Reactivity with Sulfoxides and Thioethers[END_REF] with the slopes and intercepts of all three plots being remarkably similar, about -3.0 x 10 -3 Å ppm -1 , and ~2.43 Å, respectively. That the vinylidene complex cis-RuCl2(PAN)(PPh3)(C=CHPh) had one of the longest RuPA bonds within all the complexes was rationalized in terms of decreased Ru to PA π-backbonding because of competition with such bonding in the Ru=C moiety and, consistent with this, the shortest RuPA bond is seen for the five-coordinate species RuCl2(PAN)[P(p-tolyl)3] 1 (no solid-state structure was obtained for the corresponding PPh3 complex), where maximum Ru to PA π-backbonding is likely be available. [START_REF] Ma | Reactions of phenylacetylene and p-tolylacetylene with a five-coordinate Ru II complex[END_REF] The question of why the shortest bond results in a higher chemical shift remained unanswered. Indirect evidence for such 'π-back-bonding' effects, based on solid-state NMR studies [START_REF] Olivieri | A Correlation between P-31 NMR Chemical-Shift Tensors in Solid Phosphates and P-O Bond Lengths -Implications for Proton-Transfer Processes[END_REF][START_REF] Un | Study of 31 P NMR Chemical Shift Tensors and Their Correlation to Molecular Structure[END_REF] and theoretical work [START_REF] Ling | Deciphering the NMR Fingerprints of the Disordered System with Quantum Chemical Studies[END_REF] on phosphoryl-containing compounds suggests that the degree of π-bonding between P-and O-atoms is a key factor in governing the [START_REF] Jameson | Absolute Shielding Scale for 31 P from Gas phase NMR Studies[END_REF] P chemical shift tensors and PO bond lengths. Thus, a DFT study on the RuCl2(PAN)(PR3)L complexes seemed highly appropriate to help quantify details of the observed, inverse NMR/structural plot of Figure 1, and this led to the collaborative research described in this current paper. A theoretical study by Bühl's group [START_REF] Ortuño | Computational Insight into 103 Rh Chemical Shift-Structure Correlations in Rhodium Bis(phosphine) Complexes[END_REF] that derived a correlation between experimental or computed RhP bond lengths and computed 103 Rh chemical shifts for a series of Rh phosphine-containing complexes also encouraged us to venture into this area; of note, an increasing RhP bond length corresponded to an increase in the 103 Rh chemical shift, which is opposite to the correlation between the 31 P chemical shift and the RuP distance seen in Figure 1.

COMPUTATIONAL DETAILS

All geometry optimizations were performed on the full molecules (neutral or charged) considered in vacuum with the Gaussian 09 package 17 at the PBE0-D3BJ level. [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF][START_REF] Grimme | Effect of the Damping Function in Dispersion Corrected Density Functional Theory[END_REF] Ru was represented by relativistic effective core potential (RECP) from the Stuttgart group and the associated basis set. [START_REF] Andrae | Energy-adjustable ab initio pseudopotentials for the second and third row transition elements[END_REF] The remaining atoms (H, C, N, O, Cl, Br, and P) were represented by a triple-ζ pcseg-2 basis set. [START_REF] Jensen | Unifying General and Segmented Contracted Basis Sets. Segmented Polarization Consistent Basis Sets[END_REF] The analytically calculated harmonic frequencies were used to evaluate the zero-point energies (ZPE), the enthalpies and the Gibbs energies of all species at p = 1 atm and a temperature of 298.15 K. The coordinates of all optimized structures and associated energies are given in the Supporting Information, SI.

Calculations of NMR shieldings, s, and related chemical shifts, d, (eq 1, see below for s ref ), were performed within the gauge-including atomic orbital (GIAO) framework using ADF 2016 [START_REF] Te Velde | Chemistry with ADF[END_REF] with the PBE0 functional on the full complexes in vacuum. Ru, P, and all atoms directly bonded to the metal and the P of interest were represented by Slater-type basis sets of triple-ζ quality (TZ2P), the remaining atoms being represented by a double-ζ quality Slater-type basis sets (DZP). Relativistic effects were treated by the 2-component zeroth order regular approximation (ZORA) including the spin-orbit contribution. [START_REF] Van Lenthe | Relativistic regular two-component Hamiltonians[END_REF] The results were found to improve with the inclusion of the spin-orbit contribution even if its magnitude is globally small (larger for the five-coordinate than the sixcoordinate complexes). [START_REF] Vícha | Relativistic Heavy-Neigbor-Atom Effects on NMR Shifts: Concepts and Trends across the Periodic Table[END_REF] The analysis of the electron density and scalar-relativistic natural localized orbitals was carried out with the NBO 6.0 program. [START_REF] Weinhold | Valency and Bonding: A Natural Bond Orbital Donor-Acceptor Perspective[END_REF] Calculated NMR shielding tensors were analyzed on the basis of the scalar-relativistic natural localized molecular orbitals (NLMOs). [START_REF] Bohmann | Natural chemical shielding analysis of nuclear magnetic resonance shielding tensors from gauge-including atomic orbital calculations[END_REF] To allow comparison between the species, all complexes were described with a consistent Lewis structure, involving three facial bonds: RuPA. RuN, and RuP. It was verified that this Lewis structure represents properly the computed electron density for all complexes.

δ= σ ref -σ 1-σ ref ≈ σ ref -σ (1)
Literature has shown that calculations of the 31 P NMR shieldings, and thus of the chemical shifts, are challenging, since their values cover a range of 2000 ppm; a 10% accuracy in the calculation was recently still considered as "excellent". [START_REF] Feindel | Phosphorus magnetic shielding tensors for transition-metal compounds containing phosphine, phosphido, and phosphinidene ligands: Insights from computational chemistry[END_REF] Even the shielding in the simplest PH3 is challenging to calculate. [START_REF] Stoychev | Efficient and Accurate Prediction of Nuclear Magnetic Resonance Shielding Tensors with Double-Hybrid Density Functional Theory[END_REF][START_REF] Sun | Accurate prediction of nuclear magnetic resonance shielding constants: towards the accuracy of CCSD(T) complete basis set limit[END_REF][START_REF] Teal | Benchmarking densityfunctional theory calculations of NMR shielding constants and spin-rotation constants using accurate coupled-cluster calculations[END_REF] The value of 328.

RESULTS

Structural Features of RuCl2(PAN)(PR3)L.

The coordinates and a 3D representation of the optimized structures of all complexes are provided in the SI. For all ligands except L = CO, the sixcoordinate complex of calculated lowest energy corresponds to the observed isomer, implying their thermodynamic preference; this is compatible with the synthetic method of the RuCl2(PAN)(PR3)L complexes; i.e. by addition of L to the five-coordinate RuCl2(PAN)PR3, in which the empty coordination site is trans to PA, or via ligand exchange within Ru precursors in solution. Of note, HCCPh, is converted to coordinated C=CHPh, [START_REF] Ma | Reactions of phenylacetylene and p-tolylacetylene with a five-coordinate Ru II complex[END_REF] and the NO + complex is synthesized via a different route. [START_REF] Da Silva | Studies on the Ru(II) monocationic complexes [RuCl2)(NO)(PN)(PPh3)]PF6, where PN = odiphenylphosphino-N,N-dimethylaniline, and R = Ph and p-X-C6H4 (X = OMe, Me, F)[END_REF] In contrast, the CO complex, trans-RuCl2(PAN)(PR3)(CO), is formed, by an irreversible reaction of CO to the powdered, five-coordinate trans-RuCl2(PAN)(PR3); calculations show that the cis isomer is energetically preferred by ∆H = 9.6 kcal mol -1 , (Table 1) suggesting that the trans isomer is a kinetic product. 33 The solvent used for recording the 31 P{ 1 H} chemical shift is listed with the associated reference. a)

CDCl3, 8 b) CD2Cl2, 4 c) CH2Cl2, 9 d) CDCl3, 7 and e) CDCl3. [START_REF] Ma | Binding of CO and NH3 at a five-coordinate Ru(II) centre in the solid state and in solution[END_REF] NA not available.

The isomeric preference for the thermodynamically preferred complexes can be rationalized considering the s-donating and -accepting properties of L. [START_REF] Tsipis | Building trans-philicity (trans-effect/trans-influence) ladders for octahedral complexes by using an NMR probe[END_REF][START_REF] Khoroshun | Electronic reorganization: Origin of sigma trans promotion effect[END_REF] Thus, the vinylidene group which is both a strong s-donor and a single face -acceptor prefers to be trans to one of the two chlorides than to PA (see Scheme 1 for labeling) by 11 kcal mol -1 , hence the preference for the cis isomer. Indeed, a chloride ligand has a weak trans influence and also enhances the back-donation of one Ru lone-pair into the empty p orbital of vinylidene Ca via 4e d-p interaction. [START_REF] See For Instance Caulton | The influence of -stabilized unsaturation and filled-filled repulsions in transition metal chemistry[END_REF] The larger donation (L ® Ru) and backdonation (Ru ® p(CaC=CHPh)) in the cis-over the trans-isomer are confirmed by the NBO analysis (Table S18 in the SI). Similar effects account for the preference for the cis over the trans isomer for L = NO + , N2, N2O (bonded via N as shown in the SI), and H2 by ∆G of 12.7, 10.7, 7.1, and 6.9 kcal mol -1 , respectively. Thus, all L ligands with -accepting ability [START_REF]It is not easy to estimate the relative -accepting ability of H2 relative to ligands like NO + , N2O and N2 but it is[END_REF] coordinate cis to PA with the notable exception of CO, whose complex is not under thermodynamic control.

When the L ligand is a modest s-donor and has no -accepting ability, the energy preference for the cis isomer is either small, or the trans isomer is slightly more stable. The five-coordinate RuCl2(PAN)(PR3) has the empty coordination site trans to PA as illustrated in Scheme 1. [START_REF] Mudalige | Molecular-Hydrogen, -Nitrogen and Monohydride Derivatives of the Structurally Characterized Dichloro(odiphenylphosphino-N,N-dimethylaniline) [tris(p-tolyl)phosphine]ruthenium(II) Complex[END_REF] If L interacts weakly with the Ru complex, the structural preference of the five-coordinate species determines that of the six-coordinate RuCl2(PAN)(PR3)L; this accounts for the preference of trans-RuCl2(PAN)(PR3) (H2O). The marginally stronger ligand H2S induces a slightly larger perturbation to the five-coordinate species, resulting in almost isoenergetic cis-and trans-isomers; this is illustrated by the preference for the cis-isomer by 3.4 or 1.1 kcal mol -1 when the monodentate phosphine is PPh3 or P(p-tolyl)3, respectively. Under kinetic control, the preference for the trans isomer in the case of the CO ligand results from a direct addition to the vacant site of the five-coordinate Ru complex.

There is a general fair agreement between calculated and observed bond lengths even if the calculated bond distances are generally slightly shorter than the distances given by X-ray diffraction values. (Tables S4 andS6-S9 in the SI). This is likely due to intermolecular interactions present in the crystal that are not included in the calculations. Geometry optimizations carried out for a single molecule in vacuum, with a computational method that includes the dispersion corrections have been shown to give bond distances that are too short. [START_REF] Moellmann | Influence of Crystal Packing on an Organometallic Ruthenium(IV) Complex Structure: The Right Distance for the Right Reason[END_REF] A noticeable discrepancy is obtained for L = NO + .

For this ligand, the solid-state structures of RuCl2(PAN)(PR3)(NO) + reveal that Ru-N(nitrosyl) bond distance depends on the nature of the monodentate PR3, as shown by values of 1.794(5), 1.870(5), and 1.738(3) Å, for P(p-tolyl)3, PPh3, and P(p-FC6H4)3, respectively (Table S8 in the SI). [START_REF] Da Silva | Studies on the Ru(II) monocationic complexes [RuCl2)(NO)(PN)(PPh3)]PF6, where PN = odiphenylphosphino-N,N-dimethylaniline, and R = Ph and p-X-C6H4 (X = OMe, Me, F)[END_REF] Geometry optimization carried out with PPh3 does not reproduce the observed, long RuN(nitrosyl) distance.

However, since the [START_REF] Jameson | Absolute Shielding Scale for 31 P from Gas phase NMR Studies[END_REF] PA chemical shifts (and other spectroscopic data) are similar for the three phosphine ligands, [START_REF] Da Silva | Studies on the Ru(II) monocationic complexes [RuCl2)(NO)(PN)(PPh3)]PF6, where PN = odiphenylphosphino-N,N-dimethylaniline, and R = Ph and p-X-C6H4 (X = OMe, Me, F)[END_REF] the Ru-N(nitrosyl) distance is probably similar in solution for the three NO + complexes.

Isotropic 31 PA Chemical Shifts. The 31 PA chemical shifts, d P A , of RuCl2(PAN)(PR3)L, were calculated for the cis-and trans-isomers, regardless of the structural preferences, in order to determine the impact of the coordination site of L on the chemical shift of PA, (Table 1 and Figure 2).

The 31 P chemical shifts of the monodentate phosphine were also calculated, they vary little with L and are less informative, the corresponding data are reported in the SI and are not discussed in detail (Table S3). Good agreement is obtained between experimental and calculated chemical shifts for PA, with the exception of the five-coordinate RuCl2(PAN)[P(p-tolyl)3] complex (#1) for which the observed d P A in CDCl3 is 81.5 ppm 1 and the calculated value is 113.8 ppm. Since the difference between experimental and calculated values is < 15 ppm for all other complexes, this larger discrepancy of ~ 32 ppm may not be due to the computational method. This five-coordinate complex, in the absence of an agostic interaction (not present in the solid state or in the computed optimized structure) likely binds a CDCl3 ligand and forms six-coordinate RuCl2(PAN)(PR3)(CDCl3) in solution.

Optimization of this complex gave a trans h 1 Cl3CD complex (#2) with a Ru … Cl distance of 2.79 Å and also a less stable cis isomer (#3); further, for this complex #2, the calculated d P A of 87.4 ppm is only 6 ppm away from the experimental value. This is considered a good indication that the NMR data for this five-coordinate complex, in fact, refer to a solvated species with a weakly coordinated CDCl3.

The agreement between calculated and experimental chemical shifts for PA in this set of complexes is good as illustrated by an R 2 of 0.964, excluding #1 from the correlation (Figure 2).

Including dynamic effects and explicit solvation of the molecular species, in addition to representation of the relativistic effects with 4-components calculations, would be required to achieve high accuracy as illustrated recently for trans-PtCl2(NHMe2)(PPh3) in water. [START_REF] Castro | Four-component relativistic 31 P NMR calculations for trans-platinum(II) complexes: importance of the solvent and dynamics in spectral simulations[END_REF] In fact, the lack of consideration of the diversity of possible conformations could account for the 13 ppm difference between the observed and calculated chemical shifts in the case of L = C=CHPh (#21). For probably related dynamic effects, the agreement between calculated and experimental 31 P chemical shifts was found to be significantly lower for the monodentate phosphine than for PA of the bidentate PAN ligand (Table S3 in the SI vs squares and red empty triangles correspond to observed and non-observed isomers, respectively. For labelling of complexes, see Table 1.

As expected, the calculated chemical shifts differ for the cis and trans isomers, but the difference between these values depends on L. With the weakly bonded ligands H2O, H2S, CH3SH, N2O, N2 and H2 the calculated d P A are similar for the two isomers with no clear trend for which isomer has the higher chemical shift. It is thus rather difficult to use the calculated chemical shift to assign the cis or trans nature of the isomer by comparison with the experimental shift. The situation is better with more strongly bonded ligands such as NO + , C=CHPh and CO (#18-23). This confirms the cis nature for the experimental NO + and C=CHPh complexes and trans nature for the CO complex.

Similarly, the calculations indicate clearly that the empty coordination site (or the weakly interacting CDCl3 solvent) is trans to PA ((#1-2).

DISCUSSION

Charge transfer from L to RuCl2(PAN)(PR3)

To understand how L influences the chemical shifts of PA, one needs first to analyze how L tunes the electronic structure of RuCl2(PAN)(PR3) around PA. A meaningful way to obtain this information is to consider that the six-coordinate complex is formed by the union of RuCl2(PAN)(PR3) and L, and study this interaction using the NBO method. [START_REF] Weinhold | Valency and Bonding: A Natural Bond Orbital Donor-Acceptor Perspective[END_REF] In an NBO vision, the bonding is described by donor/acceptor interactions; the key one is delocalization of the electron density on L, usually localized in a lone-pair according to the NBO analysis, into vacant antibonding s* type NBOs (antibonds) between Ru and the ligands of the fivecoordinate fragment. [START_REF] Weinhold | Valency and Bonding: A Natural Bond Orbital Donor-Acceptor Perspective[END_REF]40 This delocalization leads to charge transfer from L to the vacant s* NBOs;

due to the absence of symmetry of the complexes, charge transfer occurs also in antibonds that are cis and trans to L. As expected, the charge transfer increases with the s-donor strength of L but is always significantly less than 2e. Therefore, the partial occupancy of a s* NBO antibond leads to lengthening but not cleavage of the Ru-ligand bond. In term of molecular orbitals (which can always be expressed as linear combinations of NBOs), partial occupancy of a s* NBO is associated with a raise in energy of all vacant (canonical) molecular orbitals that have contributions from this s* NBO.

Here, focus on the RuPA bond is needed, and thus the partial occupancy of the s(RuPA) and s*(RuPA) NBOs were calculated (Table S17 in the SI). The partial occupancy of s(RuPA) NBO is around 1.89e for all L. In contrast, the partial occupancy of s*(RuPA), ns*(RuPA), is very sensitive to the nature of L. This partial occupancy, which is the smallest (0.15e) for the five-coordinate complex #1, increases in a modest manner (0.2 to 0.4e) for weak s-donors like H2O, H2S, CH3SH, and N2O, (#4-13) and becomes significant (0.5 to 0.65e) for strong s-donor L like NO + , C=CHPh, and CO (#18-23). Increasing ns*(RuPA) lengthens the RuPA bond (Figure S6 in the SI) although the correlation is modest (R 2 = 0.652) because several additional donor/acceptor s-and -interactions contribute to the global interaction of the two fragments.

Correlation between Isotropic Shielding at PA and Electronic Properties at PA. The following uses both the chemical shift (d P A ) and the shielding (s P A ) terms. As currently recognized, [START_REF] Viesser | The unexpected roles of σ and π orbitals in electron donor and acceptor group effects on the 13 C NMR chemical shifts in substituted benzenes[END_REF][START_REF] Halbert | Elucidating the Link between NMR Chemical Shifts and Electronic Structure in d 0 Olefin Metathesis Catalysts[END_REF][START_REF] Yamamoto | Orbital Analysis of Carbon-13 Chemical Shift Tensors Reveals Patterns to Distinguish Fischer and Schrock Carbenes[END_REF] the NMR shielding of PA is not expected to correlate well with the atomic charge on the active nucleus. This is the case also for the complexes in the present work (Figure S7 in the SI). In attempts at understanding the observed correlation between the d P A /s P A and the RuPA bond distance, a search was made for an electronic property that could correlate with these factors, and a remarkably good correlation (R 2 = 0.912, Figure 4) was found between the calculated shielding at PA, s P A , and the electron occupancy of the s*(RuPA) NBO, ns*(RuPA), defined in the previous section. Thus, the isotropic shielding and the RuPA bond distance both increase as ns*(RuPA) increases, i.e. as more electron density is transferred to the s*(RuPA) NBO. This naturally establishes a correlation between d P A /s P A and the RuPA bond distance, as observed experimentally. While this is gratifying phenomenologically, an understanding is developed in the following section. increasing number of studies focused on the relationship between the chemical shifts/shieldings and bonding features of the selected nucleus and the interested reader is referred to selected articles and reviews. [START_REF] Bohmann | Natural chemical shielding analysis of nuclear magnetic resonance shielding tensors from gauge-including atomic orbital calculations[END_REF][START_REF] Facelli | Chemical shift tensors: Theory and application to molecular structural problems[END_REF][START_REF] Widdifield | Understanding Chemical Shielding Tensors Using Group Theory, MO Analysis, and Modern Density-Functional Theory[END_REF][START_REF] Bühl | The DFT route to NMR chemical shifts[END_REF][START_REF] Kaupp | Interpretation of NMR Chemical Shifts[END_REF][START_REF] Falivene | Theoretical NMR spectroscopy of N-heterocyclic carbenes and their metal complexes[END_REF] These studies have been extremely useful in getting deeper insight on the bonding properties around the studied atom, with possible extrapolation to other properties such as reactivity. 49- 55 The flourish of recent studies for a high diversity of atoms is a clear indication that chemists wish to learn more NMR information. The literature is abundant for main group elements and only selected references are given. Carbon [START_REF] Gordon | Metathesis Activity Encoded in the Metallacyclobutane Carbon-13 NMR Chemical Shift Tensors[END_REF][START_REF] Marchione | C-13 NMR Spectroscopy of N-Heterocyclic Carbenes Can Selectively Probe sigma Donation in Gold(I) Complexes[END_REF][START_REF] Gordon | NMR chemical shift analysis decodes olefin oligo-and polymerization activity of d 0 group 4 metal complexes[END_REF][START_REF] Gordon | Metal alkyls programmed to generate metal alkylidenes by a-H abstraction: prognosis from NMR chemical shift[END_REF][START_REF] Gordon | π-Bond Character in Metal-Alkyl Compounds for C-H Activation: How, When, and Why?[END_REF][START_REF] Gordon | Carbon-13 NMR Chemical Shift: A Descriptor for Electronic Structure and Reactivity of Organometallic Compounds[END_REF][START_REF] Zilm | Low-temperature natural-abundance Carbon-13 NMR spectroscopy of matrix-isolated Species. The anisotropy of the shielding tensor in ethylene[END_REF][START_REF] Azofra | Quantifying electronic similarities between NHC-gold(I) complexes and their isolobal imidazolium precursors[END_REF] and silicon [START_REF] Cavalieri | Unusual Silicon-29 Chemical-Shift Anisotropies in Three-Membered Rings[END_REF] were considerably studied. More recently aluminium, [START_REF] Lam | Understanding Trends in Al-27 Chemical Shifts and Quadrupolar Coupling Constants in Chloroalkyl Aluminum [AlClx(Me)(3-x) ](n = 1 or 2) Compounds[END_REF] yttrium, 60 oxygen, 61 selenium, 62 and tellurium 63 were also considered. Phosphorus has been studied in systems with and without a transition metal. [START_REF] Duchamp | Nature of the carbon-phosphorus double Bond and the carbon-phosphorus triple bond as studied by solid-state NMR[END_REF][START_REF] Wu | Unusual 31 P Chemical Shielding Tensors in Terminal Phosphido Complexes Containing a Phosphorus-Metal Triple Bond[END_REF][START_REF] Eichele | Phosphorus-31 chemical shift tensors of phosphinidene ligands in ruthenium carbonyl clusters: a 31 P single[END_REF][67][68] For instance, the nature of the carbon-phosphorus double and triple bonds was interpreted by solid-state NMR. [START_REF] Duchamp | Nature of the carbon-phosphorus double Bond and the carbon-phosphorus triple bond as studied by solid-state NMR[END_REF] The unusually high chemical shift of a metal-bound phosphido group was shown by calculations to be associated with a high paramagnetic term. [START_REF] Wu | Unusual 31 P Chemical Shielding Tensors in Terminal Phosphido Complexes Containing a Phosphorus-Metal Triple Bond[END_REF] Solid-state NMR studies and analysis of the chemical shift tensor have been combined on a variety of phosphine-containing complexes using computational methods similar to the one used in the current work. [START_REF] Eichele | Phosphorus-31 chemical shift tensors of phosphinidene ligands in ruthenium carbonyl clusters: a 31 P single[END_REF][67][68] The key information is the chemical shift anisotropy (CSA), which can be measured by solidstate NMR and calculated by quantum methods. It is represented by a second rank tensor whose components in the principal axes frame are d11, d22, d33 (d11 ≥ d22 ≥ d33); the isotropic chemical shift is the average of these three components (equation 2) and a similar quantity in terms of shieldings (s11 ≤ s22 ≤ s33) is used below. A measure of the anisotropy is represented by the span W -difference between d11 and d33.

δ iso = 1 3 ( δ 11 +δ 22 + δ 33 ) (2)
For all complexes computed in this work, one of the principal axes of the shielding tensor is approximatively along the RuPA direction and the other two are approximatively perpendicular to it.

(Table S11, and Figures S3,S4, and S5 in the SI). Thus, although the complexes have no symmetry elements, the directions of the principal axes are characteristic of an octahedral type coordination as established for phosphine complexes of higher symmetry. 67 As a consequence of the absence of symmetry elements, the two principal axes that are perpendicular to the RuPA "axial" direction do not have specific orientation relative to the four ligands in the "equatorial plane". Their orientations change with L and its cis or trans position relative to PA. However, for essentially all complexes, the two axes perpendicular to the RuPA direction are associated with the most deshielded components, i.e.

σ 11 and σ 22 . Of note, however, for cis NO + (#19), and to a lesser extent cis C=CHPh (#21), the principal axes associated with the s22 and s33 components deviate from being perpendicular and parallel, respectively, to the RuPA direction. Fortunately, these deviations do not impact significantly on the analysis.

In a nonrelativistic formalism, the shielding can be split into diamagnetic and paramagnetic (usually deshielding) contributions. 69 Within the GIAO framework, the calculations show that the diamagnetic contribution for PA in all species of Table 1 is essentially constant and is also isotropic (Table S12 in the SI). Thus, the variation in chemical shifts as a function of L must be attributed to the paramagnetic contribution, which includes also the spin-orbit term in this work. This is nicely illustrated by the almost perfect correlation between siso and s p+so with a slope of 0.995 and R 2 of 0.9997 (Figure 5 left and Table S12 in the SI). Analysis of the terms σ 11 p+SO , σ 22 p+SO , and σ 33 p+SO , which are the paramagnetic plus spin-orbit shieldings in the principal axis frame, is thus sufficient to explain the behavior of s p+so , and thus of siso.

Calculations show that σ 11 p+SO and σ 22 p+SO vary significantly with L whereas σ 33 p+SO is more constant (Table S12, in the SI). Therefore, this latter term cannot be responsible for the variation of siso. Further, the correlations between siso and σ 11 p+SO or σ 22 p+SO , independently, are reasonable (R 2 of 0.930 and 0.884, respectively, Figure S1 in the SI). However, the slopes of these correlations are similar, suggesting equal contributions. Indeed, the correlation between siso and 1/3(σ 11

p+SO + σ 22 p+SO )
with an excellent R 2 of 0.9863 and a slope of 0.844 shows that the paramagnetic contributions of the two most deshielded terms to the averaged paramagnetic shielding s p+so quantitatively describes the changes of the isotropic shielding siso. We now analyze the reasons for this behavior. An analysis of the shielding principal components s11, s22 and s33, was performed in terms of molecular orbitals (MOs). To understand the origin of the contributions of the various MOs to the paramagnetic terms 70 a convenient approach is that of the 'orbital rotation model'. [START_REF] Zilm | Low-temperature natural-abundance Carbon-13 NMR spectroscopy of matrix-isolated Species. The anisotropy of the shielding tensor in ethylene[END_REF]71,72 Focus being on phosphorus, all MOs having contributions on the PAC and RuPA bonds are of interest. The orbital rotation model is illustrated in Scheme 2 with selected contributions of the phosphorus orbitals to these bonds. For instance, the 3py orbital contributes to MOs having PAC bonding character. Applying the magnetic field along the x direction can be thought of as rotating the 3py orbital around the x axis by 90°, this yielding a rotated orbital aligned with the z axis. Thus, the occupied MOs with 3py orbital character overlap, after rotation, with an empty MO with s*(RuPA) character. This coupling results in a paramagnetic contribution at PA associated with the x direction. An analogous development can be done with the PA 3pz orbital, which contributes to MOs having RuPA bonding character. Applying the magnetic field along the x direction can be thought of as rotating the 3pz orbital around the x axis, this yielding a rotated orbital along the y direction. The occupied MOs with 3pz character overlap, after rotation, with an empty orbital that is antisymmetrical relative to the xz plane such as the schematic orbital shown on the right-hand side of Scheme 2. The corresponding coupling between these orbitals also yields a paramagnetic term associated with the x axis. The magnitude of the paramagnetic contributions increases with increasing overlap between the occupied and empty MOs and decreasing energy gap. Similar schemes apply to the magnetic field along the y direction (rotation around y). In this case, the 3px and 3pz orbitals of phosphorus contribute to a paramagnetic term associated with the y direction. Finally, the magnetic field along the z direction is associated with contributions from MOs with 3px and 3py characters and a resulting paramagnetic field associated with the z direction. The paramagnetic terms s11 and s22 are associated with the magnetic field along the x and y axes, while s33 is associated with that along z axis. The contributions of the MOs were calculated for the three principal components. For all of them, a large number of small contributions is found, making interpretation of the results in term of bonding pattern challenging. Presenting the results in a more compact way is thus preferable. This is achieved with the localization procedure provided by the Natural Chemical Shift (NCS) analysis, implemented in the NBO method; this analysis describes the diamagnetic and paramagnetic contributions, notably in terms of lone-pairs and bonds. Comparison between complexes #1-23 is possible if a similar Lewis structure is selected for all systems.

Therefore, the NCS analysis was carried out using the Lewis structure that was also used to calculate the occupancy of the s*(RuPA) NBO. This procedure reduces considerably the number of terms; four emerge: PAC and RuPA bonds, the Ru lone-pairs and PA core orbitals that contribute roughly in comparable amounts to the paramagnetic terms of s11, s22 and s33 (Tables S13-S15). Unfortunately, none of these four contributions correlate well with the full paramagnetic term of each principal component for all complexes #1-23. However, it is sufficient to consider the contributions of the Mentioned earlier was that empty MOs with s*(RuPA) NBO character are expected to be low in energy when ns*(RuPA) is small. This is the case for #1-2 where ns*(RuPA) is < 0.25e (Table S17); this leads to large PAC bond deshielding paramagnetic contributions to s11 and s22, for #1-2, as indeed calculated (Tables S13 andS14). Following the same analysis, there is much less deshielding contribution from the PAC bonds in the case of #18-23; for these systems, ns*(RuPA) is large (around 0.5-0.6e), and in fact the calculations give a shielding contribution. A related analysis can be done for the contribution of the s(RuPA) NBO to the paramagnetic terms of s11 and s22. A schematic representation of an empty orbital, which can couple with the rotated s(RuPA) NBO, is shown on the right-hand side of Scheme 2 in the case of L = CO. Such an empty orbital is at low energy when L has -accepting character. This accounts for a significant deshielding paramagnetic contribution originating from the RuPA bond for ligands like NO + , C=CHPh and CO (#18-23). Thus, the contributions from PAC and RuPA bonds have opposite consequences on the shielding of PA but the contribution from the PAC bonds varies more across the whole set #1-23, which results in PA being the most deshielded in the complexes #1-2 and the least in the complexes #18-23. This also explains why ns*(RuPA) appears to be a good compact descriptor for the PA chemical shift. In addition, ns*(RuPA)

correlates well with the chemical shifts for the entire series #1-23, (Figure 4) which is probably due to additional, non-identified, small contributions.

The 31 P chemical shift is thus highly sensitive to the coordination sphere of the phosphine complex, but cannot be related to the donor/acceptor property of the phosphine acting as a response for a change of L. Analysis of the 31 P chemical shift is based on the occupied and empty orbital pattern of the whole complex, which is itself tuned by the coordinated ligands. However, informative is that a good correlation is shown to exist between the 31 P chemical shift and the occupancy of a s*(RuPA)

NBO, ns*(RuPA), obtained from an NBO analysis of the complex considered as the union of RuCl2(PAN)(PR3) and L. This suggests that the PA chemical shift is more strongly influenced by the s metal-ligand field than by its donor/acceptor characteristics.

Figure 1 .

 1 Figure 1. Plot of 31 PA{ 1 H} experimental data vs. experimental RuPA bond length (Å) for RuCl2(PAN) [P(p-tolyl)3] and the RuX2(PAN)(PR3)L complexes, where X = Cl or Br, and R = Ph, p-tolyl or p-FC6H4 ; L is shown in Scheme 1; various crystalline forms are included when occurring; none is RuCl2(PAN)[P(p-tolyl)3].

Figure 1 Figure 3

 13 Figure1presented the linear relationship between experimental[START_REF] Jameson | Absolute Shielding Scale for 31 P from Gas phase NMR Studies[END_REF] PA chemical shifts in solution and RuPA bond distances determined in the solid state. A similar relationship holds for computed values and even when all isomers (observed and not observed) are included, as shown in Figure3where blue squares and red triangles refer to the observed and non-observed species, respectively. However, the correlation is better when considering the observed isomers only (blue line R 2 = 0.960) than all isomers (black line R 2 = 0.857). Three complexes are far from the correlation line, like RuCl2(PAN)(P(p-tolyl)3) (#1), trans-RuCl2(PAN)(PPh3)(NO) + (#18) and also cis-RuCl2(PAN) (PPh3)(C=CHPh) (#21). Thus, the correlation between d P A and RuPA bond distances appears good, but the presence of exceptions suggests care in a more general application.

Figure 3 .

 3 Figure 3. Calculated chemical shifts d P A (ppm) vs. calculated RuPA distances (Å). The blue full

Figure 4 .

 4 Figure 4. Isotropic shielding vs. the occupancy, ns*(RuPA), of the s*(RuPA) NBO.

Figure 5 .

 5 Figure 5. Isotropic shielding siso vs. (left) s p+so and (right) (σ 11 p+SO + σ 22 p+SO )/3.

  PAC and RuPA bonds to s11 and s22 to understand why the five-coordinate species and its solvated variant (#1-2), have the most deshielded PA, and six-coordinate species with L = NO + , C=CHPh and CO, (#18-23), have the most shielded PA. In fact, the contributions of the PAC bonds to the paramagnetic terms of s11 and s22 vary considerably, and alone determine the low shieldings for #1-2 and the high ones for #18-23. An opposite, but smaller, variation is found from the contribution of the RuPA bond; thus, this contribution moderates the very large change in shielding associated with the contribution from the PAC bonds. Together, the sum of these two contributions account well qualitatively for the high and low chemical shifts of complexes #1-2 and #18-23, respectively.Nevertheless, these two terms, together, do not correlate well with the variation of the paramagnetic terms of s11 and s22 for the entire series, #1-23; additional contributions are needed. Since the contributions from the PAC and RuPA bonds are needed to define the complexes with the lowest and highest PA chemical shifts, they are worth exploring further.As suggested by the NCS analysis, and illustrated in Scheme 2, the contribution of the PAC term to the paramagnetic term along the x direction depends on the 'overlap' of the rotated PA 3py orbital with the s*(RuPA) NBO. It also depends on the difference in energy between the MOs with PAC and s*(RuPA) NBOs contributions. Qualitatively, the MOs with PAC contributions are likely to be influenced little by the coordination at the ruthenium center. In contrast, the energies of empty MOs with s*(RuPA) NBO character are likely to be significantly influenced by this coordination and use of the partial occupancy ns*(RuPA) of the s*(RuPA) NBO is proposed to establish tendencies.

Scheme 2 .

 2 Scheme 2. Schematic representation of the rotated orbital model for a magnetic field along the x direction resulting into paramagnetic couplings. (Left) 3py orbital, contributing to the PAC bond, coupled with the s*(RuPA) NBO after 90° rotation around x axis; L is not specified. (Right) 3pz orbital, contributing to the RuPA bond, coupled with a schematic vacant orbital (L = CO) after rotation by 90° around the x axis

Table 1 .

 1 Relative enthalpies (∆H) and Gibbs energies (∆G ) in kcal mol -1 for the cis and trans isomers

	of RuCl2(PAN)(PR3)L; experimental and calculated 31 PA NMR chemical shifts in ppm using the
	calculated value of 351.6 ppm for the isotropic shielding of the reference, 32 and the experimental and
	calculated RuPA distances in Å
	#	L isomer R(PR3)	∆H	∆G 31 P δexp 31 P δcalc	Ru PA expt	RuPA calc
	1	none	trans p-tolyl	-	-	81.5 a	113.8	2.170(1) 1	2.154
	2	CHCl3	trans p-tolyl -9.7 -10.5		87.4		2.194
	3	CHCl3	cis p-tolyl	0.0	0.0		68.6		2.259
	4	H2O	trans phenyl -8.0	-8.1	68.5 a	77.5 2.2305(14) 8	2.203
	5	H2O	cis phenyl	0.0	0.0		76.5		2.219
	6	CH3SH	trans phenyl	4.3	4.5		47.7		2.280
	7	CH3SH	cis phenyl	0.0	0.0	51.4 a	60.3	2.2802(8) 5	2.248
	8	H2S	trans phenyl	4.1	3.4		50.8		2.273
	9	H2S	cis phenyl	0.0	0.0	50.6 a	55.3	2.2712(6) 5	2.253
	10	H2S	trans p-tolyl	2.2	1.1		51.8		2.280
	11	H2S	cis p-tolyl	0.0	0.0	51.9 a	55.7	2.2560(4) 2,	2.252
	12	N2O	trans phenyl	8.1	7.1		68.6		2.255
	13	N2O	cis phenyl	0.0	0.0	49.1 b	60.7	NA	2.249
	14	N2	trans phenyl	8.5	10.7		58.8		2.305
	15	N2	cis phenyl	0.0	0.0	47.2 b	57.5	NA	2.263
	16	H2	trans phenyl	4.9	6.9		52.1		2.273
	17	H2	cis phenyl	0.0	0.0	49.3 a	54.8	2.2884(7) 3	2.252
	18	NO +	trans phenyl 13.0	12.7		43.5		2.460
	19	NO +	cis phenyl	0.0	0.0	30.9 c	34.3 2.3937(10) 9	2.344
	20 C=CHPh	trans phenyl 11.7	11.0		-2.5		2.515
	21 C=CHPh	cis phenyl	0.0	0.0	37.9 d	25.2	2.337(2) 7	2.304
	22	CO	trans phenyl	9.6	10.2	13.0 e	6.9	NA	2.432
	23	CO	cis phenyl	0.0	0.0		39.0		2.272

Table 1 )

 1

.

Figure

2

. Calculated vs. experimental

[START_REF] Jameson | Absolute Shielding Scale for 31 P from Gas phase NMR Studies[END_REF]

PA chemical shifts for RuCl2(PAN)(PR3)L. See Table

1

for labelling of the complexes. The reference value for calculating

[START_REF] Jameson | Absolute Shielding Scale for 31 P from Gas phase NMR Studies[END_REF]

P chemical shifts is 351.6 ppm.

[START_REF] Lantto | NMR Shielding Constants in PH3. Absolute Shielding Scale. and the Nuclear Magnetic Moment of 31 P[END_REF]

ACKNOWLEDGEMENTS

This work was initiated following a visit of OE to UBC with lectureship from the French Lecture series sponsored by the French Embassy in Canada. OE was partially supported by the Research Council of Norway through its Centre of Excellence Scheme Project (Hylleraas Centre, grant 262695) and the earlier Centre of Excellence project (CTCC (grant 179568/V30). The authors are grateful to Dr Heike Fliegl now at the Institute of Nanotechnology (INT), Karlsruhe Institute of Technology (KIT), who performed early test calculations when she was at the CTCC in Oslo and to Dr Michal Repisky from the Hylleraas centre at the Artic University of Norway for useful discussions. Presented calculations were carried out in Montpellier using the French national facilities (GenCI grant A0070810475).

AUTHOR INFORMATION

Corresponding Authors *E-mail christophe.raynaud1@umontpellier.fr *E-mail : odile.eisenstein@umontpellier.fr *E-mail: brj@chem.ubc.ca

ORCID

TABLE OF CONTENT GRAPHIC and SYNOPSIS

The experimental correlation established between RuP distances and 31 P NMR chemical shifts in RuCl2(PN)(PR3)L (L = none, H2O, H2S, CH3SH, H2, N2, N2O, NO + , C=CHPh, and CO; PN = odiphenylphosphino-N,N'-dimethylaniline) is well reproduced with DFT calculations and, through a study of the P chemical shift tensors and the paramagnetic contributions to the shielding, shown to be mostly associated with the s-donor strength of the L ligands.


Huyên Pham

Xavier Warin

Maximilien Germain

Neural networks-based backward scheme for fully nonlinear PDEs *

Keywords: Neural networks, fully nonlinear PDEs in high dimension, backward scheme. MSC Classification: 60H35, 65C20, 65M12

We propose a numerical method for solving high dimensional fully nonlinear partial differential equations (PDEs). Our algorithm estimates simultaneously by backward time induction the solution and its gradient by multi-layer neural networks, while the Hessian is approximated by automatic differentiation of the gradient at previous step. This methodology extends to the fully nonlinear case the approach recently proposed in [HPW20] for semi-linear PDEs. Numerical tests illustrate the performance and accuracy of our method on several examples in high dimension with non-linearity on the Hessian term including a linear quadratic control problem with control on the diffusion coefficient, Monge-Ampère equation and Hamilton-Jacobi-Bellman equation in portfolio optimization.

Introduction

This paper is devoted to the resolution in high dimension of fully nonlinear parabolic partial differential equations (PDEs) of the form ∂ t u + f (., ., u, D x u, D 2 x u) = 0, on [0, T ) × R d , u(T, .) = g, on R d ,

(1.1) with a non-linearity in the solution, its gradient D x u and its hessian D 2 x u via the function f (t, x, y, z, γ) defined on [0, T ] × R d × R × R d ×S d (where S d is the set of symmetric d × d matrices), and a terminal condition g.

The numerical resolution of this class of PDEs is far more difficult than the one of classical semi-linear PDEs where the nonlinear function f does not depend on γ. In fact, rather few methods are available to solve fully nonlinear equations even in moderate dimension.

• First based on the work of [START_REF] Cheridito | Second-order backward stochastic differential equations and fully nonlinear parabolic PDEs[END_REF], an effective scheme developed in [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] using some regression techniques has been shown to be convergent under some ellipticity conditions later removed by [START_REF] Tan | A splitting method for fully nonlinear degenerate parabolic PDEs[END_REF].

Due to the use of basis functions, this scheme does not permit to solve PDE in dimension greater than 5.

• A scheme based on nesting Monte Carlo has been recently proposed in [START_REF] Warin | Monte Carlo for high-dimensional degenerated Semi Linear and Full Non Linear PDEs[END_REF]. It seems to be effective in very high dimension for maturities T not too long and linearities not too important.

• A numerical algorithm to solve fully nonlinear equations has been proposed by [START_REF] Beck | Machine Learning Approximation Algorithms for High-Dimensional Fully Nonlinear Partial Differential Equations and Second-order Backward Stochastic Differential Equations[END_REF] based on the second order backward stochastic differential equations (2BSDE) representation of [START_REF] Cheridito | Second-order backward stochastic differential equations and fully nonlinear parabolic PDEs[END_REF] and global deep neural networks minimizing a terminal objective function, but no test on real fully nonlinear case is given. This extends the idea introduced in the pioneering papers [EHJ17; HJE18], which were the first serious works for using machine learning methods to solve high dimensional PDEs.

• The Deep Galerkin method proposed in [START_REF] Sirignano | DGM: A deep learning algorithm for solving partial differential equations[END_REF] based on some machine learning techniques and using some automatic differentiation of the solution seems to be effective on some cases. It has been tested in [START_REF] Al-Aradi | Solving Nonlinear and High-Dimensional Partial Differential Equations via Deep Learning[END_REF] for example on the Merton problem.

In this article, we introduce a numerical method based on machine learning techniques and backward in time iterations, which extends the proposed schemes in [START_REF] Vidales | Unbiased deep solvers for parametric PDEs[END_REF] for linear problems, and in the recent work [START_REF] Huré | Deep backward schemes for high-dimensional nonlinear PDEs[END_REF] for semi-linear PDEs. The approach in these works consists in estimating simultaneously the solution and its gradient by multi-layer neural networks by minimizing a sequence of loss functions defined in backward induction. A basic idea to extend this method to the fully nonlinear case would rely on the representation proposed in [START_REF] Cheridito | Second-order backward stochastic differential equations and fully nonlinear parabolic PDEs[END_REF]: at each time step t n of an Euler scheme, the Hessian D 2

x u at t n is approximated by a neural network minimizing some local L 2 criterion associated to a BSDE involving D x u at date t n+1 and D 2

x u. Then, the pair (u, D x u) at date t n is approximated/learned with a second minimization similarly as in the method described by [START_REF] Huré | Deep backward schemes for high-dimensional nonlinear PDEs[END_REF]. The first minimization can be implemented with different variations but numerical results show that the global scheme does not scale well with the dimension. Instability on the D 2

x u calculation rapidly propagates during the backward resolution. Besides, the methodology appears to be costly when using two optimizations at each time step. An alternative approach that we develop here, is to combine the ideas of [START_REF] Huré | Deep backward schemes for high-dimensional nonlinear PDEs[END_REF] and the splitting method in [START_REF] Beck | Deep splitting method for parabolic PDEs[END_REF] in order to derive a new deep learning scheme that requires only one local optimization during the backward resolution for learning the pair (u, D x u) and approximating D 2

x u by automatic differentiation of the gradient computed at the previous step. The outline of the paper is organized as follows. In Section 2, we briefly recall the mathematical description of the classical feedforward approximation, and then derive the proposed neural networks-based backward scheme. We test our method in Section 3 on various examples. First we illustrate our results with a PDE involving a nonlinearity of type uD 2

x u. Then, we consider a stochastic linear quadratic problem with controlled volatility where an analytic solution is available, and we test the performance and accuracy of our algorithm up to dimension 20. Next, we apply our algorithm to a Monge-Ampère equation, and finally, we provide numerical tests for the solution to fully nonlinear Hamilton-Jacobi-Bellman equation, with non-linearities of the form |D x u| 2 /D 2 x u, arising in portfolio selection problem with stochastic volatilities.

The proposed deep backward scheme

Our aim is to numerically approximate the function u : [0, T ] × R d → R, assumed to be the unique smooth solution to the fully nonlinear PDE (1.1) under suitable conditions. This will be achieved by means of neural networks approximations for u and its gradient D x u, relying on a backward scheme and training simulated data of some forward diffusion process. Approximations of PDE in high dimension by neural networks have now become quite popular, and are supported theoretically by recent results in [START_REF] Hutzenthaler | Overcoming the curse of dimensionality in the numerical approximation of semilinear parabolic partial differential equations[END_REF] and [START_REF] Darbon | Overcoming the curse of dimensionality for some Hamilton-Jacobi partial differential equations via neural network architectures[END_REF] showing their efficiency to overcome the curse of dimensionality.

Feedforward neural network to approximate functions

We denote by d 0 the dimension of the input variables, and d 1 the dimension of the output variable. A (deep) neural network is characterized by a number of layers L + 1 ∈ N \ {1, 2} with m , = 0, . . . , L, the number of neurons (units or nodes) on each layer: the first layer is the input layer with m 0 = d, the last layer is the output layer with m L = d 1 , and the L -1 layers between are called hidden layers, where we choose for simplicity the same dimension m = m, = 1, . . . , L -1.

A feedforward neural network is a function from R d0 to R d1 defined as the composition

x ∈ R d -→ A L • • A L-1 • . . . • • A 1 (x) ∈ R . (2.1)
Here A , = 1, . . . , L are affine transformations:

A 1 maps from R d0 to R m , A 2 , . . . , A L-1 map from R m to R m , and A L maps from R m to R d1 , represented by A (x) = W x + β ,
for a matrix W called weight, and a vector β called bias term, : R → R is a nonlinear function, called activation function, and applied component-wise on the outputs of A , i.e., (x 1 , . . . , x m ) = ( (x 1 ), . . . , (x m )). Standard examples of activation functions are the sigmoid, the ReLu, the Elu, tanh. All these matrices W and vectors β , = 1, . . . , L, are the parameters of the neural network, and can be identified with an element θ ∈ R Nm , where

N m = L-1 =0 m (1+m +1 ) = d 0 (1+m)+m(1+m)(L-2)+m(1+d 1 )
is the number of parameters. We denote by N d0,d1,L,m the set of all functions generated by (2.1) for θ ∈ R Nm .

Forward-backward representation

Let us introduce a forward diffusion process

X t = X 0 + t 0 µ(s, X s )ds + t 0 σ(s, X s )dW s , 0 ≤ t ≤ T, (2.2)
where µ is a function defined on [0, T ] × R d with values in R d , σ is a function defined on [0, T ] × R d with values in M d the set of d × d matrices, and W a d-dimensional Brownian motion on some probability space (Ω, F, P) equipped with a filtration F = (F t ) 0≤t≤T satisfying the usual conditions. The process X will be used for the simulation of training data in our deep learning algorithm, and we shall discuss later the choice of the drift and diffusion coefficients µ and σ, see Remark 2.3. Let us next denote by (Y, Z, Γ) the triple of F-adapted processes valued in R × R d ×S d , defined by

Y t = u(t, X t ), Z t = D x u(t, X t ), Γ t = D 2 x u(t, X t ), 0 ≤ t ≤ T.
(2.3) By Itô's formula applied to u(t, X t ), and since u is solution to (1.1), we see that (Y, Z, Γ) satisfies the backward equation:

Y t = g(X T ) - T t µ(s, X s ).Z s + 1 2 tr(σσ (s, X s )Γ s ) -f (s, X s , Y s , Z s , Γ s ) ds - T t σ (s, X s )Z s .dW s , 0 ≤ t ≤ T.
(2.4)

Remark 2.1 This BSDE does not uniquely characterize a triple (Y, Z, Γ) contrarily to the semilinear case (without a non-linearity with respect to Γ) in which proper assumptions on the equation coefficients provide existence and uniqueness for a solution couple (Y, Z). In the present case at least two options can be used to estimate the Γ component:

• Rely on the 2BSDE representation from [START_REF] Cheridito | Second-order backward stochastic differential equations and fully nonlinear parabolic PDEs[END_REF] which extends the probabilistic representation of [START_REF] Pardoux | Adapted solution of a backward stochastic differential equation[END_REF] for semilinear equations to the fully nonlinear case. It is the approach used by [START_REF] Beck | Machine Learning Approximation Algorithms for High-Dimensional Fully Nonlinear Partial Differential Equations and Second-order Backward Stochastic Differential Equations[END_REF] with a global large minimization problem, as in [START_REF] Han | Solving high-dimensional partial differential equations using deep learning[END_REF].

• Compute the second order derivative by automatic differentiation. This is the point of view we adopt in this paper together with a local approach solving several small optimization problems. In this way, we provide an extension of [START_REF] Huré | Deep backward schemes for high-dimensional nonlinear PDEs[END_REF] to cover a broader range of nonlinear PDEs.

Algorithm

We now provide a numerical approximation of the forward backward system (2.2)-(2.4), and consequently of the solution u (as well as its gradient D x u) to the PDE (1.1). We start from a time grid π = {t i , i = 0, . . . , N } of [0, T ], with t 0 = 0 < t 1 < . . . < t N = T , and time steps ∆t i := t i+1 -t i , i = 0, . . . , N -1. The time discretization of the forward process X on π is then equal (typically when µ and σ are constants) or approximated by an Euler scheme:

X ti+1 = X ti + µ(t i , X ti )∆t i + σ(t i , X ti )∆W ti , i = 0, . . . , N -1,
where we set ∆W ti := W ti+1 -W ti (by misuse of notation, we keep the same notation X for the continuous time diffusion process and its Euler scheme). The backward SDE (2.4) is approximated by the time discretized scheme

Y ti Y ti+1 -µ(t i , X ti ).Z ti + 1 2 tr σσ (t i , X ti Γ ti -f (t i , X ti , Y ti , Z ti , Γ ti ) ∆t i -σ (t i , X ti )Z ti .∆W ti ,
that is written in forward form as

Y ti+1 F (t i , X ti , Y ti , Z ti , Γ ti , ∆t i , ∆W ti ), i = 0, . . . , N -1, (2.5) with F (t, x, y, z, γ, h, ∆) := y -f (t, x, y, z, γ)h + z σ(t, x)∆, (2.6) f (t, x, y, z, γ) := f (t, x, y, z, γ) -µ(t, x).z - 1 2 tr σσ (t, x)γ .
The idea of the proposed scheme is the following. Similarly as in [START_REF] Huré | Deep backward schemes for high-dimensional nonlinear PDEs[END_REF], we approximate at each time t i , u(t i , .) and its gradient D x u(t i , .), by neural networks x ∈ R d → (U i (x; θ), Z i (x; θ)) with parameter θ that are learned optimally by backward induction: suppose that Ûi+1 := U i+1 (.; θ * i+1 ), Ẑi+1 := Z i+1 (.; θ * i+1 ) is an approximation of u(t i+1 , .) and D x u(t i+1 , .) at time t i+1 , then θ * i is computed from the minimization of the quadratic loss function:

Li (θ) = E Ûi+1 -F (t i , X ti , U i (X ti ; θ), Z i (X ti ; θ), D Ẑi+1 (T (X ti+1 )), ∆t i , ∆W ti ) 2
where T is a truncation operator such that T (X) is bounded for example by a quantile of the diffusion process and D Ẑi+1 stands for the automatic differentiation of Ẑi+1 . The idea behind the truncation is the following. During one step resolution, the estimation of the gradient is less accurate at the edge of the explored domain where samples are rarely generated. Differentiating the gradient gives a very oscillating Hessian at the edge of the domain. At the following time step resolution, these oscillations propagate to the gradient and the solution even if the domain where the oscillations occur is rarely attained. In order to avoid these oscillations, a truncation is achieved, permits to avoid that the oscillations of the neural network fit in zone where the simulations propagate scarcely to areas of importance. This truncation may be necessary to get convergence on some rather difficult cases. Of course this truncation is only valid if the real Hessian does not varies too much.

The intuition for the relevance of this scheme to the approximation of the PDE (1.1) is the following. From (2.3) and (2.5), the solution u to (1.1) should approximately satisfy

u(t i+1 , X ti+1 ) F (t i , X ti , u(t i , X ti ), D x u(t i , X ti ), D 2 x u(t i , X ti ), ∆t i , ∆W ti ).
Suppose that at time t i+1 , Ûi+1 is an estimation of u(t i+1 , .). Recalling the expression of F in (2.6), the quadratic loss function at time t i is then approximately equal to

Li (θ) E u(t i , X ti ) -U i (X ti ; θ) + D x u(t i , X ti ) -Z i (X ti ; θ) σ(t i , X ti )∆W ti -∆t i f (t i , X ti , u(t i , X ti ), D x u(t i , X ti ), D 2 x u(t i , X ti )) -f (t i , X ti , U i (X ti ; θ), Z i (X ti ; θ), D Ẑi+1 (T (X ti+1 ))) 2 .
By assuming that f has small non-linearities in its arguments (y, z, γ), say Lipschitz, possibly with a suitable choice of µ, σ, the loss function is thus approximately equal to

Li (θ) (1 + O(∆t i ))E u(t i , X ti ) -U i (X ti ; θ) 2 + O(∆t i )E D x u(t i , X ti ) -Z i (X ti ; θ) 2 + O(|∆t i | 2 ).
Therefore, by minimizing over θ this quadratic loss function, via stochastic gradient descent (SGD) based on simulations of (X ti , X ti+1 , ∆W ti ) (called training data in the machine learning language), one expects the neural networks U i and Z i to learn/approximate better and better the functions u(t i , .) and D x u(t i , ) in view of the universal approximation theorem for neural networks. The rigorous convergence of this algorithm is postponed to a future work.

To sum up, the global algorithm is given in Algo 1 in the case where g is Lipschitz and the derivative can be analytically calculated almost everywhere. If the derivative of g is not available, it can be calculated by automatic differentiation of the neural network approximation of g.

Algorithm 1 Algorithm for fully non-linear equations.

1: Use a single deep neural network (U N (.; θ), Z N (.; θ)) ∈ N d,1+d,L,m and minimize (by SGD)

     LN (θ) := E U N (X t N ; θ) -g(X t N ) 2 + ∆t N -1 d E Z N (X t N ; θ) -Dg(X t N ) 2 θ * N ∈ arg min θ∈R Nm LN (θ). 2: U N = U N (.; θ * N ), and set Z N = Z N (.; θ * N ) 3: for i = N -1, . . . , 0 do 4:
Use a single deep neural network (U i (.; θ), Z i (.; θ)) ∈ N d,1+d,L,m for the approximation of (u(t i , .), D x u(t i , .)), and compute (by SGD) the minimizer of the expected quadratic loss function

     Li (θ) := E U i+1 (X ti+1 ) -F (t i , X ti , U i (X ti ; θ), Z i (X ti ; θ), D Ẑi+1 (T (X ti+1 )), ∆t i , ∆W ti ) 2 θ * i ∈ arg min θ∈R Nm Li (θ).
(2.7)

5:

Update:

U i = U i (.; θ * i ), and set Z i = Z i (.; θ * i ).
Remark 2.2 Several alternatives can be implemented for the computation of the second order derivative. A natural candidate would consist in choosing to approximate the solution u at time t i by a neural network U i and estimate Γ i as the iterated automatic differentiation D 2 x U i . However, it is shown in [START_REF] Huré | Deep backward schemes for high-dimensional nonlinear PDEs[END_REF] that choosing only a single neural network for u and using its automatic derivative to estimate the Z component degrades the error in comparison to the choice of two neural networks U, Z. A similar behavior has been observed during our tests for this second order case and the most efficient choice was to compute the derivative of the Z network. This derivative can also be estimated at the current time step t i instead of t i+1 . However this method leads to an additional cost for the neural networks training by complicating the computation of the automatic gradients performed by Tensorflow during the backpropagation. It also leads numerically to worse results on the control estimation, as empirically observed in Table 5 and described in the related paragraph "Comparison with an implicit version of the scheme". For this reason, we decided to apply a splitting method and evaluate the Hessian at time t i+1 . For this reason, we decided to apply a splitting method and evaluate the Hessian at time t i+1 .

Remark 2.3 The diffusion process X is used for the training simulations in the stochastic gradient descent method for finding the minimizer of the quadratic loss function in (2.7), where the expectation is replaced by empirical average for numerical implementation. The choice of the drift and diffusion parameters are explained in Section 3.1.

Numerical results

We first construct an example with different non-linearities in the Hessian term and the solution. We graphically show that the solution is very well calculated in dimension d = 1 and then move to higher dimensions. We then use an example derived from a stochastic optimization problem with an analytic solution and show that we are able to accurately calculate the solution. Next, we consider the numerical resolution of the Monge-Ampère equation, and finally, give some tests for a fully nonlinear Hamilton-Jacobi-Bellman equation arising from portfolio optimization with stochastic volatilities.

Choice of the algorithm hyperparameters

We describe in this paragraph how we choose the various hyperparameters of the algorithm and explain the learning strategy.

• Parameters of the training simulations: the choice of the drift coefficient is typically related to the underlying probabilistic problem associated to the PDE (for example a stochastic control problem), and should drive the training process to regions of interest, e.g.., that are visited with large probability by the optimal state process in stochastic control. In practice, we can take a drift function µ(.) equal to the drift associated to some a priori control. This choice of control could be an optimal control for a related problem for which we know the solution, or could be the control obtained by the first iteration of the algorithm. The choice of the diffusion coefficient σ is also important: large σ induces a better exploration of the state space, but as we will see in most of examples below, it gives a scheme slowly converging to the solution with respect to the time discretization and it generates a higher variance on the results. Moreover, for the applications in stochastic control, we might explore some region that are visited with very small probabilities by the optimal state process, hence representing few interest. On the other hand, small σ means a weak exploration, and we might lack information and precision on some region of the state space: the solution calculated at each time step is far more sensitive to very local errors induced by the neural network approximation and tends to generate a bias. Therefore a trade off has to be found between rather high variance with slow convergence in time and fast convergence in time with a potential bias. We also refer to [START_REF] Nüsken | Solving high-dimensional Hamilton-Jacobi-Bellman PDEs using neural networks: perspectives from the theory of controlled diffusions and measures on path space[END_REF] for a discussion on the role of the diffusion coefficient.

In practice and for the numerical examples in the next section, we test the scheme for different σ and by varying the number of time steps, and if it converges to the same solution, one can consider that we have obtained the correct solution. We also show the impact of the choice of the diffusion coefficient σ.

• Parameters of truncation: Given the training simulations X, we choose a truncation operator T p indexed by a parameter p close to 1, so that T p (X t ) corresponds to a truncation of X t at a given quantile φ p . In the numerical tests, we shall vary p between 0.95 and 0.999.

• Parameters of the optimization algorithm over neural networks: In the whole numerical part, we use a classical Feedforward network using layers with m neurons each and a tanh activation function, the output layer uses an identity activation function. At each time step the resolution of equation (2.7) is achieved using a mini-batch with 1000 training trajectories. The training and learning rate adaptation procedure is the following:

• Every 40 inner gradient descent iterations, the loss is checked on 10000 validation trajectories.

• This optimization sequence is repeated with 200 outer iterations for the first optimization step at date t N = T and only 100 outer iterations at the dates t i with i < N .

• An average of the loss calculated on 10 successive outer iterations is performed. If the decrease of the average loss every 10 outer iterations is less than 5% then the learning rate is divided by 2.

The optimization is performed using the Adam gradient descent algorithm, see [START_REF] Kingma | A Method for Stochastic Optimization[END_REF]. Notice that the adaptation of the learning rate is not common with the Adam method but in our case it appears to be crucial to have a steady diminution of the loss of the objective function. The procedure is also described in [START_REF] Chan-Wai-Nam | Machine learning for semi linear PDEs[END_REF] and the chosen parameters are similar to this article. At the initial optimization step at time t N = T , the learning rate is taken equal to 1E -2 and at the following optimization steps, we start with a learning rate equal to 1E -3.

During time resolution, it is far more effective to initialize the solution of equations (2.7) with the solution (U, Z) at the next time step. Indeed the previously computed values at time step t i+1 are good approximations of the processes at time step t i if the PDE solution and its gradient are continuous. All experiments are achieved using Tensorflow [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF]. In the sequel, the PDE solutions on curves are calculated as the average of 10 runs. We provide the standard deviation associated to these results. We also show the influence of the number of neurons on the accuracy of the results.

A non-linearity in uD 2

x u

We consider a generator in the form

f (t, x, y, z, γ) = ytr(γ) + y 2 + 2y 2 -2y 4 e -(T -t) ,
and g(x) = tanh

d i=1 xi √ d
, so that an analytical solution is available:

u(t, x) = tanh d i=1 x i √ d e -T -t 2 .
We fix the horizon T = 1, and choose to evaluate the solution at t = 0 and x = 0.5 1 I d √ d (here 1I d denotes the vector in R d with all components equal to 1), for which u(t, x) = 0.761902 while its derivative is equal to 1.2966. This initial value x is chosen such that independently of the dimension the solution is varying around this point and not in a region where the tanh function is close to -1 or 1.

The coefficients of the forward process used to solve the equation are (here

I d is the identity d × d-matrix) σ = σ √ d I d , µ = 0,
and here the truncation operator is chosen equal to

T p (X 0,x t ) = min max[x -σ √ tφ p , X 0,x t ], x + σ √ tφ p ,
where φ p = N -1 (p), with N is the CDF of a unit centered Gaussian random variable.

In the numerical results, we take p = 0.999 and m = 20 neurons. We first begin in dimension d = 1, and show in Figure 1 As the dimension increases, we have to increase the value of σ of the forward process. In dimension 3, the value σ = 0.5 gives high standard deviation in the result obtained as shown on Figure 3, while in dimension 10, see Figure 4, we see that the value σ = 1 is too low to give good results. We also clearly notice that in 10D, a smaller time step should be used but in our test cases we decided to consider a maximum number of time steps equal to 160.

Convergence of u depending on σ

Convergence of Dxu (first component) depending on σ On this simple test case, the dimension is not a problem and very good results are obtained in dimension 20 or above with only 20 neurons and 2 layers.

A linear quadratic stochastic test case.

In this example, we consider a controlled process X = X α with dynamics in R d according to

dX t = (AX t + Bα t )dt + Dα t dW t , 0 ≤ t ≤ T, X 0 = x,
where W is a real Brownian motion, the control process α is valued in R, and the constant coefficients

A ∈ M d , B ∈ R d , D ∈ R d .
The quadratic cost functional to be minimized is

J(α) =E T 0 X t QX t + α 2 t N dt + X T P X T ,
where P , Q are non negative d × d symmetric matrices and N ∈ R is strictly positive. The Bellman equation associated to this stochastic control problem is:

∂u ∂t + inf a∈R (Ax + Ba).D x u + a 2 2 tr(DD D 2 x u) + x Qx + N a 2 = 0, (t, x) ∈ [0, T ) × R d , u(T, x) = x P x, x ∈ R d ,
which can be rewritten as a fully nonlinear equation in the form (1.1) with

f (t, x, y, z, γ) = x Qx + Ax.z - 1 2 |B z| 2 tr(DD γ) + 2N
.

An explicit solution to this PDE is given by

u(t, x) = x K(t)x,
where K(t) is non negative d × d symmetric matrix function solution to the Riccati equation:

K + A K + KA + Q - KBB K N + D KD = 0, K(T ) = P.
We take T = 1. The coefficients of the forward process used to solve the equation are

σ = σ √ d I d , µ(t, x) = Ax.
In our numerical example we take the following parameters for the optimization problem:

A = I d , B = D = 1I d , Q = P = 1 d I d , N = d
and we want to estimate the solution at x = 1I d .

In this example, the truncation operator (indexed by p between 0 and 1 and close to 1) is as follows:

T p (X x t ) = min max xe Ât -σ e 2 Ât - 1 2 Â φ p , X x t , xe Ât + σ e 2 Ât - 1 2 Â φ p ,
where φ p = N -1 (p), Â is a vector so that Âi = A ii , i = 1, ..., d, 1 is a unit vector, and the square root is taken componentwise.

On Figure 5 we give the solution of the PDE with d = 1 using σ = 1.5 obtained for two dates: at t = 0.5 and at t close to zero. We observe that we have a very good estimation of the function value and a correct one of the Γ value at date t = 0.5. The precision remains good for Γ close to t = 0 and very good for u and D x u.

Y at date t = 0.5.

Z at date t = 0.5 Γ at date t = 0.5 Y at date t = 0.006125. Z at date t = 0.006125 Γ at date t = 0.006125

Figure 5: Test case linear quadratic 1D using 160 time steps, σ = 1.5, p = 0.999, 100 neurons.

On Figure 6, we give the results obtained in dimension d = 1 by varying σ. For a value of σ = 2, the standard deviation of the result becomes far higher than with σ = 0.5 or 1.

Convergence of u depending on σ.

Standard deviation of u On Figure 7, for d = 3, we take a quite low truncation factor p = 0.95 and observe that the number of neurons to take has to be rather high. We have also checked that taking a number of hidden layers equal to 3 does not improve the results. On Figure 8, for d = 3, we give the same graphs for a higher truncation factor. As we take a higher truncation factor, the results are improved by taking a higher number of neurons (100 in the figure below). On Figure 9, we observe in dimension 7 the influence of the number of neurons on the result for a high truncation factor p = 0.999. We clearly have a bias for a number of neurons equal to 50. This bias disappears when the number of neurons increases to 100.

Convergence with 50 neurons

Convergence with 100 neurons Figure 9: Convergence in 7D of the linear quadratic case, 2 layers, p = 0.999.

On Figure 10, for d = 7, we check that influence of the truncation factor appears to be slow for higher dimensions.

Figure 10: Function value convergence in 7D of the linear quadratic case with 2 layers, 100 neurons, testing p, using σ = 2 .

Finally, we give results in dimension 10, 15 and 20 for p = 0.999 on Figures 11,12. We observe that the number a neurons with 2 hidden layers has to increase with the dimension but also that the increase is rather slow in contrast with the case of one hidden layer as theoretically shown in [START_REF] Pinkus | Approximation theory of the MLP model in neural networks[END_REF]. For σ = 5 we had to take 300 neurons to get very accurate results.

10D 15D

Figure 11: Function value convergence in 10D and 15D of the linear quadratic case with 2 layers, p = 0.999.

Figure 12: Function value convergence in 20D of the linear quadratic case with 2 layers, p = 0.999.

Monge-Ampère equation

Let us consider the parabolic Monge-Ampère equation

∂ t u + det(D 2 x u) = h(x), (t, x) ∈ [0, T ] × R d , u(T, x) = g(x), (3.1)
where det(D 2 x u) is the determinant of the Hessian matrix

D 2 x u. It is in the form (1.1) with f (t, x, γ) = det(γ) -h(x).
We test our algorithm by choosing a C 2 function g, then compute G = det(D 2

x g), and set h := G -1. Then, by construction, the function

u(t, x) = g(x) + T -t,
is solution to the Monge-Ampère equation (3.1). We choose g(x) = cos(

d i=1 x i / √ d)
, and we shall train with the forward process X = x 0 + W , where W is a d-dimensional Brownian motion. On this example, we use neural networks with 3 hidden layers, d + 10 neurons per layer, and we do not need to apply any truncation to the forward process X. Actually, we observe that adding a truncation worsens the results. For choosing the truncation level, we first test the method with no truncation before decreasing the quantile parameter p. In the Monge-Ampère case the best results are obtained without any truncation. It may be caused by the oscillation of the Hessian.

The following table gives the results in dimension d = 5, 15, and for T = 1. 1: Estimate of u(0, x 0 = 1 5 ) on the Monge Ampere problem (3.1) with N = 120. Average and standard deviation observed over 10 independent runs are reported.

Portfolio selection

We consider a portfolio selection problem formulated as follows. There are n risky assets of uncorrelated price process P = (P 1 , . . . , P n ) with dynamics

dP i t = P i t σ(V i t ) λ i (V i t )dt + dW i t , i = 1, . . . , n, where W = (W 1 , . . . , W n ) is a n-dimensional Brownian motion, b = (b 1 , . . . , b n
) is the rate of return of the assets, λ = (λ 1 , . . . , λ n ) is the risk premium of the assets, σ is a positive function (e.g. σ(v) = e v corresponding to the Scott model), and V = (V 1 , . . . , V n ) is the volatility factor modeled by an Ornstein-Uhlenbeck (O.U.) process

dV i t = κ i [θ i -V i t ]dt + ν i dB i t , i = 1, . . . , n, (3.2)
with κ i , θ i , ν i > 0, and

B = (B 1 , . . . , B n ) a n-dimensional Brownian motion, s.t. d < W i , B j > = δ ij ρ ij dt, with ρ i := ρ ii ∈ (-1, 1
). An agent can invest at any time an amount α t = (α 1 t , . . . , α n t ) in the stocks, which generates a wealth process X = X α governed by

dX t = n i=1 α i t σ(V i t ) λ i (V i t )dt + dW i t .
The objective of the agent is to maximize her expected utility from terminal wealth:

E U (X α T )] ← maximize over α
It is well-known that the solution to this problem can be characterized by the dynamic programming method (see e.g. [START_REF] Pham | Continuous-time Stochastic Control and Optimization with Financial Applications[END_REF]), which leads to the Hamilton-Jacobi-Bellman for the value function on [0,

T ) × R × R n :      ∂ t u + n i=1 κ i (θ i -v i )∂ vi u + 1 2 ν 2 i ∂ 2 vi u = 1 2 R(v) (∂ x u) 2 ∂ 2 xx u + n i=1 ρ i λ i (v i )ν i ∂ x u∂ 2 xvi u ∂ 2 xx u + 1 2 ρ 2 i ν 2 i (∂ 2 xvi u) 2 ∂ 2 xx u u(T, x, v) = U (x), x ∈ R, v ∈ R n , with a Sharpe ratio R(v) := |λ(v)| 2 , for v = (v 1 , . . . , v n ) ∈ (0, ∞) n .
The optimal portfolio strategy is then given in feedback form by α * t = â(t, X * t , V t ), where â = (â 1 , . . . , ân ) is given by âi

(t, x, v) = - 1 σ(v i ) λ i (v i ) ∂ x u ∂ 2 xx u + ρ i ν i ∂ 2 xvi u ∂ 2 xx u , (t, x, v = (v 1 , . . . , v n )) ∈ [0, T ) × R × R n ,
for i = 1, . . . , n. This Bellman equation is in the form (1.1) with

f (t, x, y, z, γ) = n i=1 κ i (θ i -v i )z i + 1 2 ν 2 i γ ii - 1 2 R(v) z 2 0 γ 00 - n i=1 ρ i λ i (v i )ν i z 0 γ 0i γ 00 + 1 2 ρ 2 i ν 2 i (γ 0i ) 2 γ 00 , for x = (x, v) ∈ R n+1 , z = (z 0 , . . . , z n ) ∈ R n+1 , γ = (γ ij ) 0≤i,j≤n ∈ S n+1
, and displays a high non-linearity in the Hessian argument γ.

The truncation operator indexed by a parameter p is chosen equal to

T p (X 0,x t ) = min max[x + µt -σ √ tφ p , X 0,x t ], x + µt + σ √ tφ p ,
where φ p = N -1 (p), N is the CDF of a unit centered Gaussian random variable. We use neural networks with 2 hidden layers and d + 10 neurons per layer. We shall test this example when the utility function U is of exponential form: U (x) = -exp(-ηx), with η > 0, and under different cases for which closed-form solutions are available:

(1) Merton problem. This corresponds to a degenerate case where the factor V , hence the volatility σ and the risk premium λ are constant, so that the generator of Bellman equation reduces to

f (t, x, y, z, γ) = - 1 2 |λ| 2 z 2 γ , (t, x, y, z) ∈ [0, T ] × R × R × R, (3.3)
with explicit solution given by:

u(t, x) = e -(T -t) |λ| 2 2 U (x), âi = λ i ησ .
We train with the forward process

X k+1 = X k + λ∆t k + ∆W k , k = 0, . . . , N, X 0 = x 0 .
(2) One risky asset: n = 1. A quasi-explicit solution is provided in [START_REF] Zariphopoulou | A solution approach to valuation with unhedgeable risks[END_REF]: We test our algorithm with λ(v) = λv, λ > 0, for which we have an explicit solution:

u(t, x, v) = U (x)w(t, v), with w(t, v) = exp - 1 2 T t R( V t,v s )ds L 1-ρ 2 where V t,v
w(t, v) = exp -φ(t) v 2 2 -ψ(t)v -χ(t) , (t, v) ∈ [0, T ] × R,
where (φ, ψ, χ) are solutions of the Riccati system of ODEs:

φ -2κφ -ν 2 (1 -ρ 2 )φ 2 + λ 2 = 0, φ(T ) = 0, ψ -(κ + ν 2 (1 -ρ 2 )φ)ψ + κθφ = 0, ψ(T ) = 0, χ + κθψ - ν 2 2 (-φ + (1 -ρ 2 )ψ 2 ) = 0, χ(T ) = 0,
with κ = κ + ρνλ, and explicitly given by (see e.g. Appendix in [START_REF] Schöbel | Stochastic volatility with an Ornstein-Uhlenbeck process and extension[END_REF])

φ(t) = λ 2 sinh(κ(T -t)) κ cosh(κ(T -t)) + κ sinh(κ(T -t)) ψ(t) = λ 2 κθ κ cosh(κ(T -t)) -1 κ cosh(κ(T -t)) + κ sinh(κ(T -t)) χ(t) = 1 2(1 -ρ 2 ) ln cosh(κ(T -t)) + κ κ sinh(κ(T -t)) - 1 2(1 -ρ 2 ) κ(T -t) -λ 2 (κθ) 2 κ2 sinh(κ(T -t)) κ cosh(κ(T -t)) + κ sinh(κ(T -t)) -(T -t) -λ 2 (κθ) 2 κ κ3 cosh(κ(T -t)) -1 κ cosh(κ(T -t)) + κ sinh(κ(T -t)) ,
with κ = κ 2 + 2ρνλκ + γ 2 λ 2 . We train with the forward process

X k+1 = X k + λθ∆t k + ∆W k , k = 0, . . . , N -1, X 0 = x 0 , V k+1 = V k + ν∆B k , k = 0, . . . , N -1, V 0 = θ.
(3) No leverage effect, i.e., ρ i = 0, i = 1, . . . , n. In this case, there is a quasi-explicit solution given by

u(t, x, v) = U (x)w(t, v), with w(t, v) = E exp - 1 2 T t R(V t,v s )ds , (t, v) ∈ [0, T ] × R n , (3.4)
where V t,v is the solution to (3.2), starting from v at time t. We test our algorithm with λ i (v) = λ i v i , λ i > 0, i = 1, . . . , n, v = (v 1 , . . . , v n ), for which we have an explicit solution given by

w(t, v) = exp - n i=1 φ i (t) v 2 i 2 + ψ i (t)v i + χ i (t) , (t, v) ∈ [0, T ] × R n , φ i (t) = λ 2 i sinh(κ i (T -t)) κ i sinh(κ i (T -t)) + κi cosh(κ i (T -t)) ψ i (t) = λ 2 i κ i θ i κi cosh(κ i (T -t)) -1 κ i sinh(κ i (T -t)) + κi cosh(κ i (T -t)) χ i (t) = 1 2 ln cosh(κ i (T -t)) + κ i κi sinh(κ i (T -t)) - 1 2 κ i (T -t) -λ 2 i (κ i θ i ) 2 κ2 i sinh(κ i (T -t)) κi cosh(κ i (T -t)) + κ i sinh(κ i (T -t)) -(T -t) -λ 2 (κ i θ i ) 2 κ i κ3 i cosh(κ i (T -t)) -1 κi cosh(κ i (T -t)) + κ i sinh(κ i (T -t))
,

with κi = κ 2 i + ν 2 i λ 2 i .
We train with the forward process

X k+1 = X k + ∆W k , k = 0, . . . , N -1, X 0 = x 0 , V i k+1 = V i k + ν i ∆B i k , k = 0, . . . , N -1, V i 0 = θ i , with < W, B i > t = 0.
Merton Problem. We take η = 0.5, λ = 0.6, T = 1, N = 120, and σ(v) = e v . We plot the neural networks approximation of u, D x u, D 2

x u, α (in blue) together with their analytic values (in orange). For comparison with Figures 6 and7, we report the error on the gradient and the initial control. In practice, after empirical tests, we choose p = 0.98 for the truncation. Table 3: Estimate of u(0, x 0 = 1, θ) on the One Asset problem with stochastic volatility (d = 2). Average and standard deviation observed over 10 independent runs are reported. The exact solution is -0.53609477.

No Leverage in Scott model. In the case with one asset (n = 1), we take η = 0.5, λ = 1.5, θ = 0.4, ν = 0.2, κ = 1. For all tests we choose T = 1, N = 120, and σ(v) = e v . In practice, after empirical tests, we choose p = 0.95 for the truncation. In the case with four assets (n = 4, d = 5), we take η = 0.5, λ = 1.5 1.1 2. 0.8 , θ = 0.1 0.2 0.3 0.4 , ν = 0.2 0.15 0.25 0.31 , κ = 1. 0.8 1.1 1.3 .

In the case with seven assets (n = 7, d = 8) we take η = 0.5, λ = 1.5 1.1 2. 0.8 0.5 1.7 0.9 , θ = 0.1 0.2 0.3 0.4 0.25 0.15 0.18 , ν = 0.2 0.15 0.25 0.31 0.4 0.35 0.22 , κ = 1. 0.8 1.1 1.3 0.95 0.99 1.02 .

In the case with nine assets (n = 9, d = 10), we take η = 0.5, λ = 1.5 1.1 2. 0.8 0.5 1.7 0.9 1. 0.9 , θ = 0.1 0.2 0.3 0.4 0.25 0.15 0.18 0.08 0.91 , ν = 0.2 0.15 0.25 0.31 0.4 0.35 0.22 0.4 0.15 , κ = 1. 0.8 1.1 1.3 0.95 0.99 1.02 1.06 1.6 .

Hamilton-Jacobi-Bellman equation from portfolio optimization is a typical example of full-nonlinearity in the second order derivative, and the above results show that our algorithm performs quite well up to dimension d = 8, but gives a high variance in dimension d = 10.

Comparison with an implicit version of the scheme. As explained in Remark 2.2, an alternative option for the estimation of the Hessian is to approximate it by the automatic differentiation of the current neural network for the Z component. It corresponds to the replacement of D Ẑi+1 (T (X ti+1 )) by DZ i (T (X ti )); θ) in (2.7). An additional change has to be made to the method for it to work. At the last optimization step (for time step t 0 = 0), we notice empirically that the variable Γ 0 is not able to properly learn the initial Hessian value at all. Therefore for this last step we use variables Y 0 , Z 0 and an explicit estimation of the second order derivative given by D Ẑ1 (T (X t1 )). We see in Table 5 that the results for the Merton problem are very similar to the ones from Table 2 for the splitting scheme but with a worse estimation of the Hessian and optimal control (the error is multiplied by around 1.5). When we tested this implicit scheme on the Monge Ampere problem we also faced computational problems during the optimization step of Tensorflow. The numerical computation of the gradient of the objective function for the backpropagation step, more precisely for the determinant part, often gives rise to matrix invertibility errors which stops the algorithm execution. For these two reasons, we focused our study on the explicit scheme.

Average Std

True value Relative error (%) u(0, x 0 = 1) -0.50572 0.00034 -0.50662 0.18 D x u(0, x 0 = 1) 0.25091 0.00067 0.25331 0.95 α(0, x 0 = 1) 0.85254 0.01956 0.80438 5.99

Table 5: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 in the Merton problem (3.3) with implicit estimation of the Hessian. Average and standard deviation (Std) observed over 10 independent runs are reported

Comparison with the 2BSDE scheme of [START_REF] Beck | Machine Learning Approximation Algorithms for High-Dimensional Fully Nonlinear Partial Differential Equations and Second-order Backward Stochastic Differential Equations[END_REF]. We conclude this paper with a comparison of our algorithm with the global scheme of [START_REF] Beck | Machine Learning Approximation Algorithms for High-Dimensional Fully Nonlinear Partial Differential Equations and Second-order Backward Stochastic Differential Equations[END_REF], called Deep 2BDSE. The tests below concern the Merton problem (3.3) but similar behavior happens on the other examples with stochastic volatilities. This scheme was implemented in the original paper only for small number of time steps (e.g. N = 30). Thus we tested this algorithm on two discretizations, respectively with N = 20 and N = 120 time steps, as shown in Figure 15, for T = 1 where we plotted the learning curve of the Deep BSDE method. These curves correspond to the values taken by the loss function during the gradient descent iterations. For this algorithm the loss function to minimize in the training of neural networks is defined as the mean L 2 error between the generated Y N value and the true terminal condition g(X N ). We observe that for this choice of maturity T = 1 the loss function oscillates during the training process and does not vanish. As a consequence the Deep 2BSDE does not converge in this case. Even when decreasing the learning rate, we noticed that we cannot obtain the convergence of the scheme.

However, the Deep 2BSDE method does converge for small maturities T , as illustrated in Table 6 with T = 0.1 and different values for the number of time steps N . Nevertheless, even if the value function is well approximated, the estimation of the gradient and control did not converge (the corresponding variance is very large), in comparison with our scheme whereas the gradient is very well approximated and the control is quite precise. We also have a much smaller variance in the results. Table 7 shows the results obtained by our method with T = 0.1 in order to compare it with the performance of [START_REF] Beck | Machine Learning Approximation Algorithms for High-Dimensional Fully Nonlinear Partial Differential Equations and Second-order Backward Stochastic Differential Equations[END_REF]. It illustrates the limitations of the global approach and justifies our introduction of a local method. Table 6: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 in the Merton problem (3.3) with maturity T = 0.1 for the [START_REF] Beck | Machine Learning Approximation Algorithms for High-Dimensional Fully Nonlinear Partial Differential Equations and Second-order Backward Stochastic Differential Equations[END_REF] scheme. Average and standard deviation observed over 10 independent runs are reported.

N Averaged value Standard deviation Theoretical value Relative error (%) u(0, x 0 = 1) 5 -0.59564 0.01136 -0.59571 0.01 u(0, x 0 = 1) 10 -0.59550 0.00037 -0.59571 0.04 u(0, x 0 = 1) 20 -0.59544 0.00054 -0.59571 0.04 D x u(0, x 0 = 1) 5 0.29848 0.00044 0.29786 0.21 D x u(0, x 0 = 1) 10 0.29842 0.00084 0.29786 0.19 D x u(0, x 0 = 1) 20 0.29785 0.00054 0.29786 0.001 α(0, x 0 = 1) 5 0.82322 0.01014 0.80438 2.34 α(0, x 0 = 1) 10 0.85284 0.07565 0.80438 6.02 α(0, x 0 = 1) 20 0.84201 0.09892 0.80438 4.68

Table 7: Estimate of the solution, its derivative and the optimal control at the initial time t = 0 in the Merton problem (3.3) with maturity T = 0.1 for our scheme. Average and standard deviation observed over 10 independent runs are reported.

Figure 2 :

 2 Figure 1: A single valuation run for test case one 1D using 160 time steps, σ = 2., p = 0.999, 20 neurons, 2 layers.

Figure 3 :Figure 4 :

 34 Figure 3: Convergence in 3D of the case one, number of neurons par layer equal to 20, 2 layers, p = 0.999.

Figure 6 :

 6 Figure 6: Convergence in 1D of the linear quadratic case, number of neurons par layer equal to 50, 2 layers, p = 0.999.

Figure 7 :

 7 Figure 7: Convergence in 3D of the linear quadratic case, 2 layers, testing the influence of the number of neurons, truncation p = 0.95.

Figure 8 :

 8 Figure 8: Convergence in 3D of the linear quadratic case, 2 layers, testing the influence of the number of neurons, truncation p = 0.99.

s

  is the solution to the modified O.U. model d Vs = κ(θ -Vs ) -ρνλ( Vs ) ds + νdB s , s ≥ t, Vt = v.

  Figure 13: Estimates of the solution and its derivatives on the Merton problem (3.3) using 120 time steps.

Figure 14 :

 14 Figure 14: Estimates of the optimal control α on the Merton problem (3.3).

Figure 15 :

 15 Figure 15: Learning curve in logarithmic scale for the scheme [BEJ19] on the Merton problem (3.3) with N = 20 times steps on the left and N = 120 time steps on the right. The maturity is T = 1 . 10000 gradient descent iterations were conducted.

Table 4 :

 4 Estimate of u(0, x 0 = 1, θ) on the No Leverage problem (3.4). Average and standard deviation observed over 10 independent runs are reported.

	Dimension d Averaged value Standard deviation Relative error (%) Theoretical solution
	2	-0.49980	0.00073	0.35	-0.501566
	5	-0.43768	0.00137	0.92	-0.441765
	8	-0.38720	0.00363	1.96	-0.394938
	10	-0.27920	0.05734	1.49	-0.275092

This work is supported by FiME, Laboratoire de Finance des Marchés de l'Energie, and the "Finance and Sustainable Devel-