Bahjat El Rez

Jiawen Liu

Virginie Béreau

Carine Duhayon

Yuki Horino

Takayoshi Suzuki

Laurent Coolen
email: coolen@insp.jussieu.fr

Jean-Pascal Sutter
email: sutter@lcc-toulouse.fr

Concomitant Emergence of Circularly Polarized Luminescence and Single-Molecule Magnet behavior in chiral-at-metal Dy complex

Circularly polarized luminescence (CPL) was evidenced for the first time in single crystals of a chiralat-metal Dy(III) single-molecule magnet. This CPL is opposite for the enantiomers and was found to appear and become stronger at lower temperatures when the relaxation of the magnetization for the Dy centers is slower. The luminescence dissymmetry factor (g lum up to +/-0.18) was found to have little dependence on the emission angle which permitted observing a similar CPL for an ensemble of micro-crystals of random orientations. The investigated Dy(III) compounds consist of bimetallic ZnLn units assembled by a Schiff base ligand, control of the stereochemistry of the Ln center was achieved using chiral [3-(trifluoromethylhydroxymethylene)-camphorate] -as anionic ligands. CPL was not found for the Eu(III) homologues.

Graphical

INTRODUCTION

Single-molecule magnets (SMM) and single-chain magnets (SCM) are the focus of a tremendous interest because of their prospective interest as active magnetic materials in future devices. These discrete and 1-D spin architectures are characterized by a possible blocking of their magnetization below a critical temperature which confers them bistability property. [START_REF] Gatteschi | Molecular Nanomagnets[END_REF] Remarkable progress has been achieved in the design of such species and in the understanding of their physics. An emerging target concerns their implementation as active materials in demonstrators, an aspect that requires addressing the readout of the orientation of their magnetization. The synergy that can exist between magnetic and optical properties offers interesting possibilities.

The interplay between light and magnetism has remained a continuous research interest since the seminal report by Faraday in 1845 [START_REF] Faraday | Experimental researches in electricity.-Nineteenth series[END_REF] of a rotation by a magnetic field of the plane of polarization of light propagating in a dielectric material. More generally, magneto-optical effects rely on a change of light polarization or intensity when propagating through a medium that is altered by the presence of a magnetic field. [START_REF] Haider | A Review of Magneto-Optic Effects and Its Application[END_REF] An interesting situation arises when the propagating medium is chiral and magnetized which results in the emergence of novel properties called magneto-chiral dichroism (MChD) in absorption and magneto-chiral birefringence in refraction. [START_REF] Wagnière | The influence of a static magnetic field on the absorption coefficient of a chiral molecule[END_REF][START_REF] Barron | Magneto-chiral birefringence and dichroism[END_REF][START_REF] Kleindienst | Interferometric detection of magnetochiral birefringence[END_REF][START_REF] Atzori | Magneto-Chiral Dichroism: A Playground for Molecular Chemists[END_REF] The simultaneous breaking of space and time-inversion symmetries leads to a slight difference in the optical absorption coefficient of the medium when the magnetic polarization of the material is parallel or antiparallel to the propagation of non-polarized light, (i.e. ↑↑ -↑↓ ≠ 0 ). [START_REF] Rikken | Pure and cascaded magnetochiral anisotropy in optical absorption[END_REF] The observed magneto-chiral dichroism will reverse sign on changing the enantiomer of the chiral medium. The same situation applies for light emitted by a luminescent chiral material. [START_REF] Rikken | Observation of magneto-chiral dichroism[END_REF][START_REF] Taniguchi | Strong magnetochiral dichroism for visible light emission in a rationally designed paramagnetic enantiopure molecule[END_REF] MChD has been evidenced in chiral magnets and chiral paramagnetic compounds, and was shown to be enhanced by the intrinsic magnetic field of the medium, especially by magnetic ordering. [START_REF] Atzori | A Chiral Prussian Blue Analogue Pushes Magneto-Chiral Dichroism Limits[END_REF][START_REF] Train | Strong magneto-chiral dichroism in enantiopure chiral ferromagnets[END_REF][START_REF] Sessoli | Strong magneto-chiral dichroism in a paramagnetic molecular helix observed by hard X-rays[END_REF] Circular polarization of luminescence induced by a magnetic field (MCPL) [START_REF] Riehl | Circularly polarized luminescence spectroscopy[END_REF] is another magnetooptical effect. It corresponds to a difference in intensities of the left-handed and right-handed circularly polarized components of the emitted light. Interestingly, chiral compounds can exhibit circularly polarized luminescence (CPL) even in absence of magnetic field. [START_REF] Zinna | Lanthanide Circularly Polarized Luminescence: Bases and Applications[END_REF][START_REF] Carr | Lanthanide complexes as chiral probes exploiting circularly polarized luminescence[END_REF] Like for MChD, MCPL intensity is stronger for ferromagnetic-like materials [START_REF] Müller | Partial circular polarization of the spectral thermal emission from ferromagnetic iron[END_REF] but the effect of intra-material magnetization on CPL remains poorly addressed. CPL can be especially strong for chiral lanthanide(III) complexes in solution. [START_REF] Zinna | Lanthanide Circularly Polarized Luminescence: Bases and Applications[END_REF] Few reports have analyzed solid-state lanthanide CPL whereas some results suggest CPL may be increased in solid-state crystals as compared to solutions. [START_REF] Zinna | Circularly polarized luminescence under near-UV excitation and structural elucidation of a Eu complex[END_REF][START_REF] Tsumatori | Circularly Polarized Light from Chiral Lanthanide(III) Complexes in Single Crystals[END_REF] Opto-electronic applications such as 3D-display, open-space secure communications or anti-counterfeiting may benefit from solidstate CPL sources.

Herein we report on circularly polarized luminescence exhibited by single crystals of a chiral-at-metal Dy(III) single-molecule magnet. This CPL is opposite for the enantiomers and was found to appear and become stronger at lower temperatures when the relaxation of the magnetization for the Dy centers is slower. Such CPL was not observed for the isostructural Eu(III) complexes with nonmagnetic ground state.

RESULTS AND DISCUSSION

The Dy(III) compounds investigated consist of bimetallic ZnLn units assembled by a Schiff base ligand (H 2 L Me2 , Scheme 1). Such ZnDy derivatives are known to exhibit SMM behavior [START_REF] Costes | Effect of Ligand Substitution around the DyIII on the SMM Properties of Dual-Luminescent Zn-Dy and Zn-Dy-Zn Complexes with Large Anisotropy Energy Barriers: A Combined Theoretical and Experimental Magnetostructural Study[END_REF][START_REF] Sun | High symmetry or low symmetry, that is the question -high performance Dy(iii) single-ion magnets by electrostatic potential design[END_REF][START_REF] Oyarzabal | Rational Electrostatic Design of Easy-Axis Magnetic Anisotropy in a ZnII-DyIII-ZnII Single-Molecule Magnet with a High Energy Barrier[END_REF][START_REF] Jia | Luminescent single-molecule magnets based on lanthanides: Design strategies, recent advances and magneto-luminescent studies[END_REF] but also to be good luminophores. [START_REF] Pasatoiu | Study of the Luminescent and Magnetic Properties of a Series of Heterodinuclear[END_REF] For an electronic transition to possess an associated CPL, it must occur in the presence of a chiral force field; in other words the luminescent metal ion must be center of chirality. This was achieved using chiral [3-(trifluoromethylhydroxymethylene)-camphorate] (camph - hereafter) as anionic ligand for the Ln center. [L Me2 Zn(Cl)Dy((+)-camph) 2 (MeOH)], Λ-Dy, was formed by reaction of DyCl 3 (H 2 O) 6 , ZnL Me2 •2H 2 O, (+)-3-(trifluoroacetyl)camphor, and piperidine in 1/1/2/2 ratio in MeOH, whereas its enantiomer ∆-Dy was obtained with (-)-3-(trifluoroacetyl)camphor. Same procedure was followed to prepare the isomorphous Eu-enantiomers Λ-Eu and ∆-Eu, and the mixed metal complex Λ-Dy/Y (see Experimental section).

Crystals formed by slow evaporation of solvent (Picture S1) crystallize in monoclinic space group P2 1 2 1 2 1 and asymmetric unit comprises a single complex. The molecular structure for ∆-Dy is depicted in Figure 1 (Figure S1 for stereo-isomer Λ-Dy, Figure S2 for Λ-Eu, and Figure S3 for the mixed metal Λ-Dy/Y derivative; selected bond distances are given in Table S2). It consists of a bimetallic complex made of the Schiff base accommodating a Zn(II) ion and a Dy(III) ion in the [N2O2] and [O2O2] compartments, respectively. The Zn(II) is in a square-base pyramid-like environment with two oxygen and two nitrogen atoms of the Schiff base in the basal position and one Cl atom in the apical position. The Dy center, in addition to the four oxygen atoms of the Schiff base, is associated with two Camph -ligands in chelate interaction with the Ln ion, and one MeOH molecule. An intramolecular hydrogen bond is established between the OH group of the alcohol and the Cl atom. The polyhedral shape of the Dy coordination sphere can be described either by a distorted capped square antiprism or a distorted muffin geometry (Table S3). [START_REF] Ruiz-Martínez | Polyhedral Structures with an Odd Number of Vertices: Nine-Coordinate Metal Compounds[END_REF] The Ln is a stereogenic center and its formation is stereospecific; the Λ conformation is obtained by coordination of (+)-camph -whereas enantiomer ∆ is achieved with (-)-camph -ligands (Figure 1). The enantiomorphism of the two derivatives is confirmed by their solid state circular dichroism (CD) spectra showing symmetrical signals but of opposite sign (figure 2). The CD transitions can be assigned to π-π* intra-ligand transitions of the Schiff-base moiety (see Figure S4). These complexes were found to exhibit slow relaxation of their magnetization. AC susceptibility revealed the emergence of an out-of-phase component, χ M '', below 12 K in zero field with however, clear contribution of quantum tunneling of the magnetization (QTM) for lower temperatures. This QTM is only partly quenched by a static field (Figure S5) therefore the dilution technique was applied to investigate the relaxation characteristics. The behavior for Λ-Dy diluted (9 %) in an isostructural diamagnetic matrix of [L Me2 Zn(Cl)Y((+)-camph) 2 (MeOH)], i.e. Λ-Dy/Y, is given in Figures 3 andS6. AC susceptibility data were collected in a static field of 1 kOe found to be the optimal field to reach slowest relaxation. Analysis of the Cole-Cole plots (Figure S6d) revealed a rather narrow distribution width of the relaxation in the higher T range (10.5-7 K) with α = 0.15. For lower T the value of α steadily increased to 0.4 for 2, K which is indicative of a multi-process relaxation.

The temperature dependence of the relaxation time (Figure S6c) was obtained by analyzing the χ M '' = f(ν) for different temperatures with an extended Debye model. The linear variation of τ = f(1/T) between 8.5 and 4 K was fitted assuming a thermally activated process yielding U eff /k B = 25 K with τ 0 = 7.4x10 -6 s. A modelling over the entire temperature range (2 to 8.5 K) was possible assuming contributions from the Orbach, Raman and direct relaxation processes ( = / + 1/ ( ) + 1/( ). The values of the adjustment parameters are U eff / k B = 23 (± 3) K, τ 0 = 1 (± 2) x10 -5 s, R = 0.1 (± 0.6) s -1 , n = 4 (± 7), D = 2 (± 2) x10 4 K -1 T -4 s -1 ). The value for τ 0 is obviously too large for an Orbach process but calculation uncertainty is large because of the number of variables. Adjustment to the experimental data failed without the Orbach contribution. These results confirm the emergence of slow relaxation of the magnetization below 12 K for probing frequency of 1 kHz. The luminescence properties of individual crystals of Λ-Dy and ∆-Dy were investigated by low temperature (5 -30 K) micro-photoluminescence. We excited an area of 30 µm diameter so that single monocrystals could be observed (Picture S2). Figure 4 shows the emission spectra for a Λ-Dy crystal along left-handed and right-handed circular polarizations. Two groups of emission peaks appear at 573 and 663 nm and are attributed respectively to the 4 F 9/2 -6 H 13/2 and 4 F 9/2 -6 H 11/2 transitions of Dy(III). [START_REF] Petoud | Eu, Tb, and Dy Chiral Lanthanide Complexes with Strong Circularly Polarized Luminescence[END_REF][START_REF] Babu | Optical Properties and White-Light Emission in Dy3+-Doped Transparent Oxyfluoride Glass and Glass Ceramics Containing CaF2 Nanocrystals[END_REF] The emission intensity decreases as a function of temperature, barely any emission being observed above 30 K. For both transitions, the left-handed circularly-polarized emission is slightly stronger than the right-handed polarization. This observation was made for all observed Λ-Dy crystals. On the other hand, for all ∆-Dy crystals the right-handed emission was consistently stronger than the left-handed emission (Figure S7). The degree of circular polarization is usually quantified by the luminescence dissymmetry factor g lum (Equation 1). 15

g lum = ! " # $ % ( ! & # ) (1)
This quantity ranges between -2 and +2 and vanishes when the emission is not circularly polarized. It is typically of the order of 10 -3 -10 -2 for chiral organic molecules but can reach higher values for lanthanide complexes, [START_REF] Zinna | Lanthanide Circularly Polarized Luminescence: Bases and Applications[END_REF] with up to g lum ≈ 1.4 reported for Eu(III). [START_REF] Lunkley | Extraordinary Circularly Polarized Luminescence Activity Exhibited by Cesium Tetrakis(3-heptafluoro-butylryl-(+)-camphorato) Eu(III) Complexes in EtOH and CHCl3 Solutions[END_REF][START_REF] Kumar | Chiral supramolecular polymerization leading to eye differentiable circular polarization in luminescence[END_REF] For solutions of Dy(III) chiral complexes, generally lower g lum values of 0.01 [START_REF] Petoud | Eu, Tb, and Dy Chiral Lanthanide Complexes with Strong Circularly Polarized Luminescence[END_REF] and up to ≈ 0.4 [START_REF] Dickins | Synthesis, Time-Resolved Luminescence, NMR Spectroscopy, Circular Dichroism and Circularly Polarised Luminescence Studies of Enantiopure Macrocyclic Lanthanide Tetraamide Complexes[END_REF] have been achieved. From our spectra, using equation 1, we extract the dissymmetry factor as a function of the wavelength (fig. 4(b)) for the two samples at 5 K. Over the two emission spectral intervals, the dissymmetry factor shows unambiguous opposite signs for the two samples, demonstrating the effect of the ligand chirality on the CPL. The dissymmetry factor is quite lower in absolute value for the ∆-Dy sample, possibly because of some measurement uncertainty. We then average the dissymmetry factor over each of the two intervals and plotted its average value as a function of the temperature in figure 4(c) and S7. In spite of some uncertainty on g lum values and some variation between the individual crystals, a general trend is observed: the degree of circular polarization decreases with temperature and vanishes around 25 K. These observations were repeated for several crystals, with a distribution of g lum values ranging between 0.04 and 0.18 for Λ-Dy and between -0.04 and -0.16 for ∆-Dy at 5 K.

It can be noticed that the appearance of CPL and the increase of g lum with decrease of temperature coincide with the emergence of slow relaxation of magnetization for these complexes. As mentioned above, CPL can be induced by the chiral structure of the luminophore but it can also be generated by a magnetic field whether applied or intrinsic to the material. For a SMM, the later results from the easy-axis orientation the individual magnetization. In order to provide additional information on the possible role of each of these sources, studies have been carried out on isomorphous Eu(III) complexes (see SI for crystal data). At low temperature, Eu(III) is essentially non-magnetic, which makes it possible to probe the contribution of structural origin. Moreover, Eu(III) is the lanthanide element for which the strongest CPL has been reported in the literature. [START_REF] Zinna | Lanthanide Circularly Polarized Luminescence: Bases and Applications[END_REF][START_REF] Petoud | Eu, Tb, and Dy Chiral Lanthanide Complexes with Strong Circularly Polarized Luminescence[END_REF][START_REF] Dickins | Synthesis, Time-Resolved Luminescence, NMR Spectroscopy, Circular Dichroism and Circularly Polarised Luminescence Studies of Enantiopure Macrocyclic Lanthanide Tetraamide Complexes[END_REF] Characteristic luminescence bands where observed below 250 K for Λ-Eu and ∆-Eu (Figure S8). It can be noticed that the most prominent emission is arising from the 5 D 0 -> 7 F 2 transition with two main lines (613 and 619 nm) for each enantiomer, whereas the intensity of the 5 D 0 -> 7 F 0 transition at 585-595 nm is weak. Their relative intensity ratio is R = 19. This feature is characteristic for a low symmetry field at the Eu(III) site, in agreement with point symmetry Cs deduced from the crystal structure (Table S3). 32,[START_REF] Judd | Hypersensitive Transitions in Rare-Earth Ions[END_REF] Down to 5 K the emission is exactly the same for the two circular polarizations. The degree of circular polarization was calculated for each peak at each temperature and g lum was always close to 0, fluctuating between -0.04 and +0.04. It can be concluded that CPL of structural origin is very low for these compounds.

Finally, we used Fourier-plane imaging to obtain the emission angular distribution in the two circular polarizations (Figure 5) for a Λ-Dy and a ∆-Dy crystal. The overall emission is distributed rather isotropically for the first crystal and has more preferred directions for the second crystal, probably due to the different shape of the crystal boundaries, but for both crystals the angular distribution is the same for the two emission lines. In spite of some significant noise, the curves clearly show that the degree of circular polarization has little dependence on the emission angle. This is confirmed by observing an ensemble of crystals of random orientations. A powder formed by crushing a single crystal and illuminated over a 30 µm diameter area showed dissymmetry factors of the same order as for a single crystal, also with isotropic angular distribution (Figure S9). There is thus little difference between the circular-polarization of a single crystal and an ensemble of randomlyoriented crystals. Beyond the obvious advantage for processing such a material, the use of a powder also solves the problem of the amount of material needed for reliable measurement.

CONCLUDING REMARKS

Concomitant rise of circularly polarized luminescence and SMM-type behavior was evidenced in the solid state for an enantiomeric pair of chiral Dy complexes. This observation strongly suggests that magnetic anisotropy of the Dy combines with structural chirality to enhance the chiral force field acting on the luminescent Ln center. The synergy between the emission and the magnetic features is supported by the lack of CPL for the Eu derivatives. This allows excluding a purely structural effect and substantiates the essential role of the magnetization blocking in the CPL observed for the Dy complexes.

It is noticeable that the observed luminescence dissymmetry factor is large with |g lum | up to 0.18. Interestingly, the same CPL was found for singles crystals and the polycrystalline powders. This is a particularly important observation as it suggests that CPL can be exploited on a crystalline powder and is not limited to a single crystal.

EXPERIMENTAL

General procedures: Reagents and solvents were obtained commercially and used without further purification. The Schiff base H 2 L Me2 and its Zn complex [ZnL Me2 ]•2H 2 O were prepared as previously reported [START_REF] Costes | Complexation of a Schiff base ligand having two coordination sites (N2O2 and O2O2) with lanthanide ions (Ln = La, Pr): an NMR study[END_REF] and recrystallized from MeOH solutions by slow evaporation of the solvents. Purity of the samples was verified by [START_REF] Gatteschi | Molecular Nanomagnets[END_REF] H NMR and elemental analysis. The two enantiomers (+) and (-) of the trifluoroacetyl-3-camphor were synthesized following the literature procedure. [START_REF] Goering | Direct determination of enantiomeric compositions with optically active nuclear magnetic resonance lanthanide shift reagents[END_REF] Purity of the samples was checked by 1 H, [START_REF] Sessoli | Strong magneto-chiral dichroism in a paramagnetic molecular helix observed by hard X-rays[END_REF] C, [START_REF] Tsumatori | Circularly Polarized Light from Chiral Lanthanide(III) Complexes in Single Crystals[END_REF] F NMR, and by specific optical rotation measurements. Elemental C, H and N analyses were performed on a Perkin-Elmer 2400 II analyzer on freshly prepared and isolated samples. IR spectra were recorded in the 4000-700 cm -1 range with a Perkin-Elmer Spectrum 100 FTIR using the ATR mode. Optical rotation measurements were performed on a Perkin-Elmer 241 polarimeter with a sodium lamp, at room temperature. [START_REF] Gatteschi | Molecular Nanomagnets[END_REF] H NMR spectra were recorded using a DPX300 Brüker spectrometer with working frequencies at 300 MHz for 1 H. Chemical shifts were referenced to the residual proton resonance of the deuterated solvents.

Magnetic studies: Magnetic measurements were carried out with a Quantum Design MPMS 5S SQUID magnetometer. Magnetic susceptibility, χ M , between 2 and 300 K was obtained with an applied field of 1000 Oe. Isothermal magnetization measurements were performed up to 5 T at 2 K. The ac susceptibility responses were recorded using an oscillating field of H AC = 3 Oe in zero field and under applied static field of 0.1T. The measurements were performed on crushed crystals from freshly isolated samples mixed to grease and held in gelatin capsules. The molar magnetic susceptibilities were corrected for sample holder and for the diamagnetic contribution of all the atoms by using Pascal's tables. 36

Solid-state Circular Dichroism (CD) measurement:

The solid-state CD spectra (collected in transmittance mode) were measured on a JASCO J-720 spectropolarimeter using KBr disks of samples. A disk was prepared from several crystals of a batch of either ∆-Dy or Λ-Dy enantiomer with crystalline KBr, according to the following procedure: A mixture of sample and KBr was well ground in an agate mortar then vacuum pressed to form a uniform disk. The disk was placed in a sample holder, which was then positioned in the cell holder of the spectropolarimeter. [START_REF] Matsushima | Transitionmetal(ii) complexes with a tripodal hexadentate ligand, 1,1,1-tris[2-aza-3-(imidazol-4-yl)prop-2enyl]ethane, exhibiting incomplete total or absolute spontaneous resolution[END_REF] For each disk, the CD signals were measured from both sides, and the average was taken. A blank disk was used to measure background spectra. The data were corrected accordingly.

Crystallographic Studies. Data collections were performed at low temperature on a Bruker Apex2 (Mo microfocus source) or an Oxford-Diffraction Gemini (Cu source) diffractometer. The structures were solved by SUPERFLIP 38 or SIR92, [START_REF] Altomare | SIR92 -a program for automatic solution of crystal structures by direct methods[END_REF] and refined by means of least-squares procedures on F or F 2 using the programs of the PC version of CRYSTALS. 40 For mixed metal compound Λ-Dy/Y best agreement was obtained for ratio 0.16 Dy/0.84. Atomic scattering factors were taken from the International tables for X-ray Crystallography. Hydrogen atoms were refined using a riding model. The absolute configuration of the compounds was determined unambiguously by Flack parameter refinement. Data are gathered in table S1. Cif have been deposited at CCDC under references 2007104 (for ∆-Dy), 2007105 (for Λ-Dy), 2018680 (for Λ-Eu), and 2033959 (for Λ-Dy/Y).

CPL studies: Crystals of size 300-800 µm were attached by a silver lacquer to the cold finger of a helium-flow cryostat (Oxford, HiRes II) and observed by a homemade microscope (0.7 numerical aperture objective, x60 magnification). The sample was illuminated through the objective by the light from a mercury lamp filtered over a 330-480 nm pass band. The fluorescence emission was collected by the same objective, filtered by a 500 nm long-pass filter and analyzed by an imaging spectrometer (Horiba Jobin-Yvon, Triax-190). Circular polarization analysis was performed by introducing a polarizer before the spectrometer entrance and comparing the +45° and -45° positions of a quarterwave plate before the polarizer.

The polarization effect of the different setup elements was corrected by introducing compensating elements. Fluorescent polymer spheres (Thermo Fisher, 200nm, 580-605 nm emission) were used for setup polarization calibration. By circular polarization analysis, an average degree of circular polarization of 0 was measured, with fluctuation in the range of ±0.02. When introducing a polarizer below the objective to polarize the emission, a 0.99 degree of linear polarization was obtained. We thus conclude that the emission polarization was preserved by the setup with 1-3 % precision.

Syntheses: Magnetic behavior is independent of the chirality therefore a full set of magnetic behavior has been collected only for Λ-Dy. The magnetic purity of the ∆ enantiomer was confirmed by susceptibility values at 300 K (found χ M T = 14.10 cm 3 mol -1 K) and AC versus T trace with frequency of 1 kHz (shown below) that was found identical to Λ-Dy.

[L Me2

DC studies:

The temperature dependence of χ M T for Λ-Dy (Fig. S5a) is in agreement with the behavior expected for an isolated Dy(III) ion. The value found at 300 K is 14.15 cm 3 mol -1 K and reaches 11.41 cm 3 mol -1 K at 4 K. Magnetization versus field at 2 K (b) reaches as value of 5.35 µ B under 5 T.

(a) (b)

AC studies: AC magnetic susceptibility as a function of temperature recorded without an applied filed is shown in (c). An out-of-phase component, χ M '', was found below 15 K and lower temperature behavior indicated other relaxation processes become active. Application of an external DC field allowed to drastically reduce these contributions but not to fully cancel them even when a field of 3 kOe (d). Since the χ M '' = f(T) for different AC frequencies obtained with H DC = 1 kOe display well defined maxima, a full data set was collected (e). The temperature dependence of the relaxation time between 3.8 and 10 K (f) was obtained by analyzing the χ M '' = f(ν) for different temperatures with an extended Debye model. The linear variation of τ =f(1/T) between 10 and 6.5 K was fitted assuming a thermally activated process yielding U eff /k B = 42 K with τ 0 = 1.4x10 -6 s (dotted blue line). A modelling over the entire temperature range

Scheme 1 :Figure 1 .

 11 Scheme 1: Sketch of the ligand system in [L Me2 Zn(Cl)Ln(camph) 2 (MeOH)]

26

 26

Figure 2 .

 2 Figure 2. Solid state CD spectra for ∆-Dy (in blue) and Λ-Dy (in red)

Figure 3 .

 3 Figure 3. AC susceptibility behavior for Λ-Dy diluted in diamagnetic [L Me2 Zn(Cl)Y((+)-camph) 2 (MeOH)].

Figure 4 .

 4 Figure 4. Temperature dependent emission of Λ-Dy crystal. (a) Purple lines: left-handed circular (LHC) polarization; green lines: right-handed circular (RHL) polarization for the Λ-Dy sample at different temperatures ; (b) CPL dissymmetry factor g lum as a function of the wavelength for the two samples at 5 K; (c) Temperature-dependent variation of the luminescence dissymmetry factor, g lum for the two samples and the two spectral bands.

Figure 5 .

 5 Figure 5. Angular distribution of the emission bands at 573-nm and 663-nm bands for a (top) Λ-Dy and (bottom) a ∆-Dy crystal at 5 K. Angular range scanned is +/-44°. The grey portions correspond to the angles not collected by the objective (outside the numerical aperture)

  Zn(Cl)Dy((+)-camph) 2 (MeOH)], Λ-Dy, and [L Me2 Zn(Cl)Dy((-)-camph) 2 (MeOH)], ∆-Dy: [ZnL Me2 ]•2H 2 O (233 mg; 0.5 mmol) and DyCl 3 •6H 2 O (188 mg; 0.5 mmol) were solubilized in 10 mL of methanol. Either the (+)-3-(trifluoroacetyl)camphor or (-)-3-(trifluoroacetyl)camphor (250 mg; 1 mmol) was added to this solution followed by piperidine (85 mg; 1 mmol). The resulting pale yellow solution was stirred at room temperature for 30 min. and then transferred in a crystallization tube for slow evaporation of the solution. Pale yellow crystals suitable for X-Ray diffraction appeared after a few days (Picture S1). They were collected and washed with methanol. Data for Λ-Dy: Yield 410 mg (71 %). IR (ATR): ν = 3251 (w), 2955 (m), 1652 (s), 1639 (s), 1535 (m), 1473 (m), 1424 (m), 1392 (w), 1293 (m), 1267 (m), 1243 (w), 1222 (vs), 1200 (vs), 1178 (s), 1126 (s), 1069 (s), 1053 (m), 1023 (m), 972 (m), 918 (w), 850 (m), 803 (w), 734 (s), 642 (m) cm -1

Figure S2 .

 S2 Figure S2. Λ-Eu : ORTEP plot (with ellipsoids at the 50% level of probability) of the asymmetric unit with numbering scheme.

Figure S5 .

 S5 Figure S5. Magnetic behaviors for Dy complex.

  (3.8 to 10 K, red line) was possible considering contributions from the Orbach and Raman relaxation processes ( The values of the adjustment parameters are U eff / k B = 53 ± 9 K, τ 0 = 0.8 ± 2 x10 -8 s, R = 0.2 ± 0.1 s -1 , n = 4.9± 0.4.

FigureFigure S6 .

 S6 FigureS5gshows the χ M '' = f(T) trace for ∆-Dy under H DC = 1 kOe and for ν = 1 kHz which confirms a behavior identical to that of the Λ-Dy enantiomer.(g)

Figure S8 .

 S8 Figure S8. Temperature dependent emission and g lum for a Λ-Eu (a, b), and a ∆-Eu (c, d) crystal (blue lines: left-handed circular (LHC) polarization; red lines: right-handed circular (RHL) polarization); (e) emission spectrum at 5 K (∆-Eu).

Figure S9 .

 S9 Figure S9. Angular distribution of the emission (bands at 573-nm and 663-nm) recorded at 5K, and temperature dependent variation of the dissymmetry factor g lum for a crystalline powder of (left) Λ-Dy and (right) a ∆-Dy crystal. The grey portions correspond to the angles not collected by the objective (outside the numerical aperture).

Figure S10 .

 S10 Figure S10.X-Ray Powder Diffraction (XRPD) at 295 K for Λ-Dy , ∆-Dy, Λ-Eu, and ∆-Eu.

Table S2 .

 S2 Selected bond distances (Å) and angles (°) for ∆and Λ-Dy, Λ-Eu, and Λ-Dy/Y complexes.

			Λ-Dy	∆-Dy	Λ-Eu	Λ-Dy/Y
	Zn1	N2	2.071(1)	2.0674(9)	2.068(2)	2.067(1)
		O2	2.0813(8)	2.0773(8)	2.081(1)	2.080(1)
		N1	2.086(1)	2.084(1)	2.0875(2)	2.086(1)
		O3	2.0907(8)	2.0904(6)	2.091(1)	2.091(1)
		Cl1	2.3116(3)	2.3091(3)	2.3077(5)	2.3087(4)
	Dy1	O5	2.323(1)	2.3214(8)	2.352(1)	2.308(1)
		O7	2.3362(9)	2.3314(8)	2.378(1)	2.324(1)
		O6	2.3556(8)	2.3529(8)	2.395(1)	2.350(1)
		O2	2.4069(8)	2.4038(8)	2.438(1)	2.395(1)
		O8	2.407(1)	2.4063(8)	2.447(1)	2.356(1)
		O3	2.4243(8)	2.4231(6)	2.4565(1)	2.4114(9
		O9	2.440(1)	2.4375(8)	2.468(1)	2.420(1)
		O4	2.5370(9)	2.5330(8)	2.550 (1)	2.524(1)
		O1	2.5583(9)	2.5536(8)	2.568(1)	2.548(1)
	C1	O1	1.435(2)	1.433(2)	1.435(2)	1.434(2)
	C2	O1	1.374(2)	1.373(1)	1.377(3)	1.381(2)
	C7	O2	1.326(1)	1.324(1)	1.327(2)	1.328(2)
	C14	O3	1.325(1)	1.322(1)	1.324(2)	1.325(2)
	C18	O4	1.375(2)	1.375(1)	1.377(2)	1.380(2)
	C19	O4	1.438(2)	1.436(1)	1.439(2)	1.442(2)
	C23	O5	1.280(2)	1.280(1)	1.277(2)	1.281(2)
	C25	O6	1.242(1)	1.238(1)	1.239(2)	1.239(2)
	C35	O7	1.273(2)	1.274(1)	1.274(2)	1.276(2)
	C37	O8	1.246(2)	1.245(1)	1.244(2)	1.245(2)
	C46	O9	1.431(2)	1.428(1)	1.432(2)	1.429(2)
	Cl1	O9	3.022(1)	3.0219(8)	3.025(1)	3.022(1)
	Cl1	H9	2.2273(3)	2.1898(3)	2.1764(4)	2.1600(4)
	O2	Zn1	80.49(3)	67.77(3)	80.95(5)	80.22(4)
		O3	80.49(3)	67.83(3)	67.20(4)	67.98(3)
	O2	Dy1
		O3
	O2

Table S3 :

 S3 SHAPE Analysis for the ML9 coordination sphere of the Ln ions.[START_REF] Gatteschi | Molecular Nanomagnets[END_REF] Program for the stereochemical analysis of molecular fragments by means of continuous shape measures and associated tools 2.1 ed.; University of Barcelona: Barcelona, 2013.

	CSAPR-9	Spherical capped square antiprism
	MFF-9	Muffin
			[ML9]	CSAPR-9 (C 4v )	MFF-9 (Cs)
			∆-Dy	1.396	1.391
			Λ-Dy	1.399	1.398
			Λ-Dy/Y	1.356	1.357
			Λ-Eu	1.545	1.457
		Ideal CSAPR	Actual polyhedron (∆-Dy)	Ideal MFF
	Figure S4. Solid-state UV-vis spectra for ZnL	Me2 and [L
								Me2 ZnL Λ-Dy ∆-Dy
		I. (a. u.)
		250	300	350	400	450	500	550
					Wavelength (nm)

Me2

Zn

(Cl)

Dy((+/-)-camph) 2 (MeOH)] complexes.

1 Llunell, M.; Casanova, D.; Cirera, J.; Alemany, P.; Alvarez, S.; SHAPE:

Acknowledgements

This work was supported by the French Research Agency (Agence Nationale de la Recherche, reference ANR-09-BLAN-0054-01). Authors are grateful to M. J.-F. Meunier and M. L. Réchignat (LCC) for technical assistance in magnetic data collection.

Conflict of Interest

The authors declare no competing financial interest.

SUPPORTING INFORMATION

CONTENTS Picture 1.

Crystals of [L

Me2

Zn (Cl)Dy((-)-camph) 2 (MeOH)], ∆-Dy. Table S1.

Crystal data for Λ-Dy and ∆-Dy.

Figure S1.

Asymmetric units for Λ-Dy and ∆-Dy with numbering scheme. Figure S2.

Λ-Eu : ORTEP plot (with ellipsoids at the 50% level of probability) of the asymmetric unit with numbering scheme.

Figure S3

Λ-Dy/Y : ORTEP plot (with ellipsoids at the 50% level of probability) of the asymmetric unit with numbering scheme.

Table S2.

Selected bond distances (Å) and angles (°) for Eu and Dy complexes Table S3.

SHAPE Analysis for the ML9 coordination sphere of the Ln ions. Crystals (size 300-800 µm) and powders (crushed crystal) on the cold stage for micro-photoluminescence studies. Figure S7.

Temperature dependent emission and g lum of a ∆-Dy crystal.

Figure S8.

Temperature dependent emission and g lum for a Λ-Eu and a ∆-Eu crystal.

Figure S9.

Angular distribution of the emission (bands at 573-nm and 663-nm) recorded at 5K, and temperature dependent variation of degree of circular polarization for a crystalline powder of Λ-Dy and a ∆-Dy crystal.


K Anselme

P Davidson

A M Popa

M Giazzon

M Liley

Lydie Ploux

Response to comment on "The interaction of cells and bacteria with surfaces structures at the nanoscale"

Response to comment on « The interaction of cells and bacteria with surfaces structures at the nanoscale."

K. Anselme 1* , P. Davidson 1 , A.M. Popa 2,3 , M. Giazzon 3 , M. Liley 3 , L. Ploux In agreement with Dr. Schilardi and collaborators, the authors are also convinced that bacteria are far from simple colloidal particles and that they are able to sense topographical features of the surface and to react accordingly, which nevertheless has not been strictly demonstrated so far. However, it cannot be expected that bacteria are completely protected from processes affecting colloidal particles, especially from dewetting processes [START_REF] Beckendorff | Quantitative and morphological analysis of biofilm formation on self-assembled monolayers[END_REF]. The author's conviction is that the experiments conducted by Dr. Schilardi and coll. in their recent article [START_REF] Díaz | Organization of Pseudomonas fluorescens on Chemically Different Nano/Microstructured Surfaces[END_REF] were affected by dewetting-related experimental artifacts. In particular, bacteria localization may have been modified by air or nitrogen drying, even in a relative humidity greater than 70% or when the sample was kept in air during a short time before immersing it again (conditions used by Schilardi and coll. for AFM imaging and for the washing procedure prior to epifluorescence imaging, respectively). A dewetting process occurs whatever the humidity of the surrounding air environment [START_REF] Elender | Wetting and dewetting of Si/SiO 2-wafers by free and lipidmonolayer covered aqueous solutions under controlled humidity[END_REF][START_REF] Sun | Effect of contact angle and humidity on evaporation of inkjet-printed colloidal drops[END_REF]. It is due to the liquid evaporation, leads to meniscus formation, convective flow, rupturing, dewetting, and capillary-bridge formation. As soon as environmental conditions allow evaporation of the colloidal suspension, colloidal particles are driven by the threephase contact lines, which will lead to their organization [START_REF] Celio | Patterned Assembly of Colloidal Particles by Confined Dewetting Lithography[END_REF][START_REF] Sun | Two-Dimensional Self-Assembly of Latex Particles in Wetting Films on Patterned Polymer Surfaces[END_REF]. Bacteria are unlikely to be able to resist to these driving forces. Our review (figure 9) [START_REF] Anselme | The interaction of cells and bacteria with surfaces structured at the nanometre scale[END_REF] displays an example of such colloidal-like behavior of bacteria. Another example, leading also to bacteria organization and alignment with grooves, is given in Figure 1. Bacteria orientated much more with grooves when colonized substrate was dried before to be immersed again in liquid. Despite the longer time the substrate stayed in air in our experiment compared to those of Dr. Schilardi and coll., this demonstrates that dewetting processes can affect bacteria organization. Finally, it can be asserted that experimental dewetting-related artifacts play a role in the results of studies addressing bacteria organization to surface, as soon as an interface between the substrate and air is created before observation. Therefore, the conclusion made by Dr. Schilardi and her group in the article [START_REF] Díaz | Organization of Pseudomonas fluorescens on Chemically Different Nano/Microstructured Surfaces[END_REF] should be rather considered as a hypothesis that still needs to be demonstrated.

Figure 1 :

 1 Figure1: Adherent Escherichia coli K12 on microstructured substrates (see[START_REF] Ploux | Opposite Responses of Cells and Bacteria to Micro/Nanopatterned Surfaces Prepared by Pulsed Plasma Polymerization and UV-Irradiation[END_REF] for more details about


N Krah

J M Létang

S Rit

Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation

Protons undergo many small angle deflections when traversing a medium, such as a patient. This effect, known as multiple Coulomb scattering (MCS), leads to degraded image resolution in proton radiography and computed tomography (CT) and to lateral spreading of the dose distribution in proton therapy. To optimally account for MCS in proton imaging, the most likely path (MLP) of a proton is estimated based on its position and propagation angle measured in front of and behind the object. In this work, we propose a functional which quantifies the likelihood of a proton trajectory and study how it can be used to model proton trajectories in a homogeneous medium. We focus on two aspects: first, we present an analytical method to quickly generate proton trajectories in a homogeneous medium based on the likelihood functional and validate it through Monte Carlo simulations. It could be used for fast generation of proton CT images without a full Monte Carlo simulation, or potentially to complement the components in a treatment planning Monte Carlo which simulate MCS. Second, by maximising the likelihood functional, we derive an expression for the MLP which is equivalent to the conventional ones reported in the literature yet computationally more convenient. Moreover, we show that the MLP is strictly a polynomial function if the protons' energy loss in the medium is approximated as a polynomial and that the orders of both are linked. We validate our MLP through Monte Carlo simulations and compare proton CT images reconstructed with our expression and with the conventional one. We find that an MLP polynomial of orders larger than five do not lead to increased spatial resolution compared to lower order expressions.

INTRODUCTION

INTRODUCTION

Protons in the therapeutic energy range (< 300 MeV) undergo a large number of small angle deflections when traversing a medium due to multiple Coulomb scattering (MCS) [START_REF] Gottschalk | Multiple Coulomb scattering of 160 MeV protons[END_REF]. In proton therapy, this leads to a successive lateral spreading of a pencil beam as it penetrates the patient tissue and its impact on the dose distribution needs to be accurately modelled by treatment planning systems. In proton radiography and computed tomography (CT), MCS is responsible for the relatively low spatial resolution compared to X-ray imaging. In these imaging techniques, protons traverse a patient at sufficiently high energy to be captured on the downstream side by a suitable detector device [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF][START_REF] Parodi | Heavy ion radiography and tomography[END_REF]. The most sophisticated systems measure each individual proton's position and propagation angle in front of and behind the patient with tracking detectors [START_REF] Schulte | Conceptual design of a proton computed tomography system for applications in proton radiation therapy[END_REF][START_REF] Scaringella | A proton Computed Tomography based medical imaging system[END_REF][START_REF] Civinini | Recent results on the development of a proton computed tomography system[END_REF][START_REF] Penfold | Geometrical optimization of a particle tracking system for proton computed tomography[END_REF][START_REF] Taylor | An experimental demonstration of a new type of proton computed tomography using a novel silicon tracking detector[END_REF]. Using this information as constraint, the most likely path (MLP) which the proton has followed across the medium can be estimated [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF][START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF][START_REF] Collins-Fekete | Developing a phenomenological model of the proton trajectory within a heterogeneous medium required for proton imaging[END_REF][START_REF] Collins-Fekete | A theoretical framework to predict the most likely ion path in particle imaging[END_REF]. The tomographic reconstruction is then performed using the MLP [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF][START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojectionthen-filtering approach[END_REF][START_REF] Penfold | Techniques in Iterative Proton CT Image Reconstruction[END_REF] and the uncertainty of the MLP estimation compared to the real trajectory is what limits the spatial resolution of the images. The conventional MLP estimation methods estimate the most likely transverse position through which a proton has passed in a certain depth and construct the MLP as a concatenation of these transverse positions by repeating the estimation in a series of depths across the medium. They exploit the fact that all MCS events are independent of each other and that the angular dispersion due to MCS is approximately Gaussian. The cited methods have been thoroughly validated with Monte Carlo simulations and one variant or another of them are commonly used in proton CT reconstruction. However, they do not provide an explicit way to quantify the likelihood of a given trajectory. This idea is elaborated from a mathematical point of view by [START_REF] Erdelyi | A comprehensive study of the most likely path formalism for protoncomputed tomography[END_REF] who discusses the concept of random fields in this context. He describes a theoretical framework to derive an MLP estimate from an action functional, however provides few explicit calculations with a direct practical application.

In this work, we introduce a likelihood functional which quantifies the likelihood of a proton trajectory based on its weighted curvature integrated along depth and investigate how it can be used to model proton trajectories in a homogeneous medium. In particular, the likelihood functional allows us to achieve two things: first, to analytically generate random proton trajectories as they would be produced by a Monte Carlo simulation in a homogeneous medium, and second, to derive an expression for the MLP by maximising the likelihood functional, as suggested by [START_REF] Erdelyi | A comprehensive study of the most likely path formalism for protoncomputed tomography[END_REF], which is equivalent to the conventional expressions. The functional disregards nuclear interactions and only takes into account MCS within the Gaussian approximation. Generating proton trajectories analytically would be of interest for quickly producing proton CT projections without running a full Monte Carlo simulation, e.g. for testing a reconstruction algorithm, or to partially replace the proton transport in a Monte Carlo code, e.g. in the context of treatment planning. We show that the MLP automatically becomes a polynomial as long as the protons' energy loss is approximated by a polynomial (which is common practice in proton CT reconstruction) and thereby create a link to the works of [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF] and [START_REF] Collins-Fekete | Developing a phenomenological model of the proton trajectory within a heterogeneous medium required for proton imaging[END_REF]. We validate our MLP expression with Monte Carlo generated data and test it in the proton CT reconstruction of a standard image quality phantom. The advantage of our expression compared to the conventional ones is that it requires fewer floating point operations to evaluate while providing the flexibility to select the desired degree of estimation accuracy by choosing the polynomial order.

MATERIAL AND METHODS

Framework of the formalism

The purpose of this work is to analytically model proton trajectories in homogeneous media. When traversing such a medium, protons undergo a large number of small angle deviations due to MCS. The likelihood of each such deflection depends on the radiation length X 0 of the material and on the proton energy E. We measure the depth coordinate u along the axis obtained by extending the proton's entrance vector across the object and denote with t the transverse position of a proton, i.e. projected onto a plane perpendicular to the u-axis. We define as trajectory the position t(u) as a function of u. The propagation angle θ of a proton projected onto this plane is given relative to the u-axis, so that θ = 0 means parallel to the u-axis. The three-dimensional trajectory is given by two transverse positions t x and t y in the two planes perpendicular to the u-axis.

We do not attempt to describe heterogeneities transverse to the beam path. Furthermore, the developed formalism accounts for multiple Coulomb scattering (MCS), but not for discrete nuclear scattering events. The main context within which we performed this investigation is proton computed tomography where nuclear events are usually filtered out and in any case not accounted for by most tomographic reconstruction algorithms.

Likelihood functional of proton trajectories

Within the Gaussian approximation of MCS, the likelihood of a change in propagation angle ∆θ = θ 2 -θ 1 over a depth interval ∆u is given by

L(∆θ) = exp - (∆θ/∆u) 2 2T ∆u = exp - (∆θ) 2 2T ∆u , ( 1
)
where T is the scattering power which depends on the proton's energy and the radiation length of the medium [START_REF] Lynch | Approximations to multiple Coulomb scattering[END_REF][START_REF] Gottschalk | Multiple Coulomb scattering of 160 MeV protons[END_REF]. The likelihood L[t] of a trajectory t can be thought of as the joint probability of a series of many deflections ∆θ i which have caused the proton to follow t(u),

L[t] = i L(∆θ i ) = i exp - (∆θ i /∆u) 2 2T i ∆u = exp - 1 2 i (∆θ i /∆u) 2 T i ∆u . (2)
Because the propagation angle θ relative to the beam axis is small, it can be taken as the derivative to the trajectory with respect to depth, θ = ∂t/∂u. Making the steps i in equation 2 above infinitesimally small and considering infinitely many of them, the sum becomes an integral and the likelihood is given as

L[t] = exp - 1 2 uout 0 (∂ 2 t(u)/∂u 2 ) 2 T (u) du = exp - 1 2 uout 0 (t (u)) 2 T (u) du , (3)
where we have replaced the finite difference ∆θ i /∆u = ((∂t/∂u) i -(∂t/∂u) i-1 )/∆u by the second derivative and used t as short hand notation. Without loss of generality, we chose to place the beginning of the trajectory at u = 0; u out denotes the depth at which the proton exits the volume. According to equation 3, the likelihood of a trajectory is given by its integrated curvature weighted by the depth-dependent scattering power T (u) to account for the proton's energy and the material's radiation length. We parametrise the scattering power as T (u) = E 2 0 /β 2 (u)p 2 (u)X 0 (u), with E 0 = 13.6 MeV/c a constant, β = v/c the proton's velocity relative to the speed of light, p the proton's momentum, and X 0 the radiation length. This is an approximate form of the expression proposed by [START_REF] Highland | Some practical remarks on multiple scattering[END_REF]. More accurately, it would contain an additional logarithmic term which is not straight forward to incorporate in the integral expression of the likelihood functional (equation 3). Some discussion on this can be found in [START_REF] Erdelyi | A comprehensive study of the most likely path formalism for protoncomputed tomography[END_REF]. For the scope of this work, we restrict ourselves to the approximate expression and leave it to future contributions to explore other parametrisations of T (u). Furthermore, we only consider homogeneous media in this work, although parametrising T (u) appropriately would allow modelling longitudinal material heterogeneities, i.e. along the u-axis [START_REF] Collins-Fekete | Extension of the FermiâĂŞEyges most-likely path in heterogeneous medium with prior knowledge information[END_REF].

Analytical generation of proton trajectories as polynomials

MCS is the most prominent physics process in a Monte Carlo simulation of protons in the therapeutic range (< 300 MeV) along with electromagnetic energy loss. To generate a proton trajectory, a Monte Carlo code has to perform many steps, each time sampling the angular deviation the proton suffers. We therefore investigated how the likelihood function (equation 3) could instead be used to analytically generate random proton trajectories as a whole. Because the likelihood function only models MCS, such an analytical generation does not account for discrete nuclear scattering events. We consider mono energetic protons propagating through a box of water. For simplicity, we will first assume that all protons have the same initial position and angle, although this is not a requirement.

We approximate proton trajectories by polynomial functions,

t(u) = N n=0 a n u n = t in + θ in u + N n=2 a n u n , (4)
so that a trajectory t is represented by an N + 1 coefficient vector a = (a 0 , a 1 , . . . a N ) T . On the right hand side, we have used the fact that

t(u = u in ) = t in = a 0 and t (u = u in ) = θ in = a 1 .
Within the framework of polynomial functions, sampling trajectories is equivalent to sampling coefficient vectors a. The second derivative t in the likelihood function (equation 3) can be calculated explicitly and is

t (u) = N n=2 n(n -1)a n u n-2 . ( 5
)
We further expand the function 1/T (u) into a polynomial,

1 T (u) ≈ L l=0 c l u l . (6)
It can be seen from figure 1 that already L = 1 gives a good approximation. Using this polynomial expression, one obtains for the weighted trajectory curvature

(t (u)) 2 T (u) = N,N,L n=2,n =2,l=0 c l n(n -1)n (n -1)a n a n u l+n+n -4 . (7)
Inserting this into the likelihood function (equation 3) and integrating the polynomial analytically yields

L poly [t] = exp   - 1 2 N,N n=2,n =2 a n a n Φ nn   (8)
with

Φ nn = L l=0 c l n(n -1)n (n -1) l + n + n -3 u l+n+n -4 out . (9)
In other words, the likelihood distribution of trajectories is approximated by an (N-1)dimensional Gaussian distribution of the coefficient vectors a. In practice, to generate a random sample of coefficients a, the following procedure is convenient: first diagonalise Φ (which is always possible because the matrix is symmetric and positive definite) to obtain Ψ = Q T ΦQ, where the columns of Q are the unit length eigenvectors of Φ and

Q T Q = 1. Defining W = diag(1/ √ λ 1 , 1/ √ λ 2 , ...1/ √ λ N ), where {λ i , i = 1, . . . N } are the eigenvalues of Φ, one has W Q T ΦQW = 1.
One can therefore sample vectors r from an (N-1)-dimensional normal distribution (unity variance) and obtain the coefficient vectors a by transformation as a = QW r.

We note that the covariance matrix of the coefficients a n is the inverse of the matrix Φ, i.e. a n a n = Φ -1 nn . This is linked to the variation of the transverse position t(u) in some depth u among an ensemble of analytically generated trajectories. Assuming for simplicity a pin-like beam impinging perpendicularly (θ in = 0) onto a box of water at t in = 0, one has t(u) = N n=2 a n u n . The average transverse position over many trajectory is zero, t(u) = 0, so that the variation of t(u) is

t 2 (u) = N n=2 a n u n N n=2 a n u n = N,N n=2,n =2 a n a n u n+n = N,N n=2,n =2 Φ -1 nn u n+n . ( 10
)
It is important to recall that the previously derived expressions correspond to only one of the two scattering planes transverse to the beam axis. A three-dimensional proton trajectory is characterised by a pair of coefficient vectors a x and a y , one for each scattering plane and each randomly sampled from the distribution described above.

By construction, the first two coefficients represent the protons' entry position and angle: t(u = 0) = t in = a 0 and t (u) = θ in = a 1 . It is worth noting that the likelihood does not depend on them (summation starts at n = 2 in equation 8) because in a homogeneous medium, the position where a proton enters the volume and its direction does not impact the likelihood of MCS. On the other hand, the coefficients a 0 and a 1 can be drawn from some distribution to generate the desired beam shape. For example, setting a x,1 = a y,1 = 0 and drawing a x,0 and a y,0 from Gaussian distributions G(a x,0

) = exp(-(a x,0 -t centre,x ) 2 /2σ 2 x )/ √ 2πσ x and G(a y,0 ) = exp(-(a y,0 -t centre,y ) 2 /2σ 2 y )/
√ 2πσ y , with standard deviation σ x = σ y , produces a symmetric parallel Gaussian pencil beam of width σ 2 beam = σ 2 x + σ 2 y centred at (t centre,x , t centre,y ).

Derivation of the most likely path

In the context of proton CT, an important concept is the most likely path, i.e. the trajectory which a proton registered by the proton CT scanner has most likely followed through the imaged object. The most sophisticated imaging systems measure positions and propagation angles of individual protons in front of and behind the object [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF]. We denote these as (t in , θ in ) and (t out , θ out ), respectively. In the conventional mathematical framework [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF], the MLP is constructed by estimating the most likely position t MLP independently at several depths u, given the measured entry and exit data. The trajectory which connects the so-obtained series of points t MLP (u) is the MLP. In the presented work, the MLP is obtained by maximising the likelihood functional L (equation 3), or by minimising L = -log L, using the Euler-Lagrange formalism commonly used in physics, e.g. to derive equations of motion.

Observing that L only depends on t , but not on t or t , we find for the infinitesimal variation of the Lagrangian (after two partial integrations)

δL = uout 0 d 2 du 2 ∂F ∂t δt du + δt ∂F ∂t uout 0 -δt d du ∂F ∂t uout 0 ∀δt ∈ C 2 , ( 11
)
where δt and δt is an arbitrary (continuous and differentiable) test function and its derivative with respect to u, respectively. We have further defined F = (t ) 2 /2T for the sake of simpler notation. Since the proton's position and angle at u = 0 and u = u out are known, δt and δt vanish there, and the last two terms in equation 11 above are zero. The minimisation requirement that δL be equal to zero for any test function δt implies

d 2 du 2 ∂F ∂t t =t MLP = 0 ⇔ d du ∂F ∂t t =t MLP = d du t MLP T = c 1 ⇔ t MLP T = c 1 u + c 0 , (12)
where c 0 and c 1 are integration constants. An expression for the MLP can be derived from the right hand side expression in equation 12 above by integration,

θ MLP (u) = t MLP (u) = θ in + c 0 u 0 dv T (v) + c 1 u 0 dv T (v)v (13) t MLP (u) = t in + θ in u + c 0 u 0 dv v 0 dw T (w) + c 1 u 0 dv v 0 dw T (w)w, (14)
where v and w are integration variables. The constants c 0 and c 1 can be determined by imposing the boundary conditions θ MLP (u = u out ) = θ out and t MLP (u = u out ) = t out . One finds

c 0 = 1 AD -BC (-B (t out -t in -θ in u out ) + D (θ out -θ in )) c 1 = 1 AD -BC (A (t out -t in -θ in u out ) -C (θ out -θ in )) , (15)
with

A = uout 0 dv T (v); C = uout 0 dv v 0 dw T (w); B = uout 0 dv T (v)v; D = uout 0 dv v 0 dw T (w)w (16)
It is interesting to note that starting from our MLP formulae (equations 13 and 14), the conventional expressions reported by [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF] and [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] can be derived. We refer to appendix A for more details.

Approximation: Describing proton energy loss by a polynomial function

In list mode proton CT reconstruction, an MLP needs to be estimated for every single recorded proton [START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF][START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF]. In its conventional representation (equation 23), several multiplications and additions of several 2 × 2 matrices need to be performed at every depth u on the reconstruction grid. The Σ and R matrices themselves can be precomputed and stored. The computational aspect was what motivated other authors to develop and implement an empirical cubic polynomial MLP [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF][START_REF] Collins-Fekete | Developing a phenomenological model of the proton trajectory within a heterogeneous medium required for proton imaging[END_REF]. The advantage is that relatively fewer multiplication operations need to be performed in that case.

In the following, we show how this can be done rigorously given that the term β 2 p 2 in T (u) is often approximated as a polynomial in terms of depth u [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. The polynomial coefficients are easily obtained by a fit to T (u) extracted from a Monte Carlo simulation of protons impinging onto a box of water. Following this strategy, we write T (u) approximately as

T (u) ≈ M m=0 b m u m . ( 17
)
An example of this approximation for different orders M is shown in figure 2 for 180 MeV and 200 MeV protons in water. The fit was performed in a depth range from 0 cm to 20 cm. It is worth observing that the goodness of the fit depends on the depth range of interest relative to the proton range because T (u) diverges too much at the end for the polynomial to represent it properly. For example, for an object of 20 cm water equivalent thickness, the discrepancy between polynomial and true scattering power is more important for 180 MeV protons (lower panel) than for 200 MeV protons (upper panel), especially beyond 15 cm. Inserting equation 17 into equation 14, the integrals can be solved analytically and the MLP

MATERIAL AND METHODS

becomes a pure polynomial expression of order M + 3:

t MLP = t in + θ in u + c 0 M m=0 b m (m + 1)(m + 2) u m+2 + c 1 M m=0 b m (m + 2)(m + 3) u m+3 . ( 18
)
The constant pre-factors c 0 and c 1 are calculated according to equation 15 using the following polynomial expressions obtained by inserting equation 17 into equation 16:

A = M m=0 b m m + 1 u m+1 out ; C = M m=0 b m (m + 1)(m + 2) u m+2 out ; B = M m=0 b m m + 2 u m+2 out ; D = M m=0 b m (m + 2)(m + 3) u m+3 out . ( 19
)
Determining the MLP thus consists in calculating the coefficients of a polynomial of order M + 3 and evaluating it in all depths across the object.

Special case: assuming constant proton energy

In a homogeneous medium, the scattering power T generally increases with depth because the protons continuously slow down. Neglecting the energy loss, one could replace T (u) by a constant. In terms of the polynomial expression for T (u) introduced in section 2.5, only one coefficient b 0 is non-zero then and

t MLP = t in + θ in u + c 0 b 0 2 u 2 + c 1 b 0 6 u 3 (20)
is of strictly cubic order. Such an expression has also been derived by [START_REF] Erdelyi | A comprehensive study of the most likely path formalism for protoncomputed tomography[END_REF] and used by [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF]. Clearly, this simplified MLP estimate still matches the measured proton position and propagation angle in front of the object because t MLP (u = 0) = t in and θ MLP (u = 0) = t MLP (u = 0) = θ in . By construction of the pre-factors c 0 and c 1 (see section 2.4), the same holds true on the exit side, i.e. t MLP (u = u out ) = t out and θ MLP (u = u out ) = t MLP (u = u out ) = θ out .

Computational efficiency of MLP estimation

How computationally intensive the polynomial expression in equation 18 is to evaluate depends on the number of arithmetic operations needed to calculate the most likely transverse position t MLP in every depth u required by a reconstruction algorithm. In this light, it is convenient to combine common powers in the two sums in equation 18 into a single sum. After an index shift (n = m + 3), one obtains

t MLP = t in + θ in u + c 0 b 0 2 u 2 + M -1 m=0 c 0 b m+1 + c 1 b m (m + 2)(m + 3) u m+3 + c 1 b M (M + 2)(M + 3) u M +3 = N =M +3 n=0 d n u n = d 0 + u (d 1 + u (d 2 + . . . u (d N -1 + d N u) . . . )) , (21)
with

d 0 = t in ; d 1 = θ in ; d 2 = c 0 b 0 2 d n,(3≤n≤M +2) = c 0 b n-2 + c 1 b n-3 n(n -1) ; d N =M +3 = c 1 b M (M + 2)(M + 3) . ( 22
)
The polynomial coefficients b m can be precomputed once before the reconstruction and the prefactors c 0 and c 1 (equation 15) need to be calculated only once per proton. Consequently, also the coefficients d n need to be calculated only once per proton. Evaluating the expanded polynomial on the right hand side of equation 21 above amounts to M + 3 multiplications and M + 3 additions, i.e. 2M + 6 floating point operations per proton and depth u. For example, only 6 or 8 floating point operations are required per depth u if the MLP is approximated by a polynomial of order M = 0 or M = 1, respectively. On the other hand, evaluating the conventional MLP expression of [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] (equation 23) for M = 5 results in at least 47 floating point operations [START_REF] Mcallister | Efficient Proton Computed Tomography Image Reconstruction using General Purpose Graphics Processing Units[END_REF], even if sparsity of the matrices R 0 and R 1 is exploited.

Verification with Monte Carlo simulations

We performed Monte Carlo simulations to verify the developed analytical methods presented in this paper using the Geant4/Gate toolkit [START_REF] Agostinelli | GEANT4 -A simulation toolkit[END_REF][START_REF] Sarrut | A review of the use and potential of the GATE Monte Carlo simulation code for radiation therapy and dosimetry applications[END_REF]. A pin-like beam of mono energetic protons was simulated impinging onto a homogeneous box of water of 50 × 50 cm 2 lateral dimension and 20 cm depth. The maximum step size was set to 1 mm and each proton's three dimensional position was stored at every step. The emstandard physics list was used to consider electromagnetic interactions only and disregard nuclear scattering as this is not included in the analytical models. We first quantified how accurate it is in general to approximate a proton trajectory resulting from MCS by a polynomial function of depth. To this end, we fitted polynomials of different orders to each simulated proton trajectory and determined the RMSE between the so approximated transverse position and the true one.

Furthermore, we quantified how closely the analytically generated trajectories would resemble those simulated with Monte Carlo by two means: first, we compared the distribution of transverse proton positions at some given depth and the beam spread as a function of depth. Second, we fitted polynomials to all simulated trajectories (see above) and compared the distribution of the fit coefficients with the analytically generated coefficient distributions (see section 2.3).

To compare the accuracy of our alternative MLP formulation with the conventional one [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] for different polynomial degrees (see section 2.5), an MLP was calculated for each proton history based on the entry and exit position and angle. The root mean square error (RMSE) between the MLP estimate and the true transverse position as a function of depth was calculated over an ensemble of 10 4 protons.

Tomographic reconstruction using the alternative MLP formulae

We reconstructed proton CT images of the CTP528 high resolution module of the Catphan phantom (The Phantom Laboratory, Salem, NY) to test our alternative MLP formulae presented in this work. Input data were simulated with Geant4/Gate using a mono energetic 200 MeV point source at 100 cm distance from iso-centre. Each proton's position and propagation angle were recorded in two tracker planes upstream and downstream from the phantom, respectively. The downstream tracker also recorded the proton's residual energy as would a calorimeter in a real proton CT scanner. The QGSP-BIC physics list was used and 3-sigma cuts applied to the exit energy and angular distribution to filter out protons which have undergone nuclear scattering [START_REF] Schulte | Density resolution of proton computed tomography[END_REF] to reflect how experimental data are typically processed. The source was rotated around the iso centre in 0.5 degree steps leading to 720 projections. A total of 2.88 × 10 8 primary protons were simulated corresponding to 225 protons • mm -2 • projection -1 . Tomographic images with 0.1 × 0.1mm 2 pixel size were reconstructed using a distance-driven binning Feldkamp-Davis-Kress algorithm [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF], where the MLP intervenes during the binning step.

RESULTS

Figure 3: Left: Proton trajectories simulated with Monte Carlo (dots) and order 6 polynomial fits to each of them (solid lines). Right: RMSE averaged over depth and over the ensemble of 10 5 proton trajectories.

RESULTS

Analytical sampling of proton trajectories

Figure 3 left shows an example of 20 proton trajectories simulated at 200 MeV in water as described in section 2.8. The dots represent the transverse positions registered by the Monte Carlo simulation. The solid lines are polynomial fits, in this case of order N = 6. The lowest two coefficients are fixed to a 0 = t in = 0 and a 1 = θ in = 0 in accordance with the proton source geometry. The right panel shows the average RMSE between the simulated trajectories and the polynomial fits (see section 2.8). Empirically, we found that the RMSE values vary approximately according to a Poisson distribution among the ensemble of 10 5 trajectories so that the uncertainties to be with the data points have the same magnitude as the RMSE values themselves.

The left panels in figure 4 show sets of 50 example trajectories analytically generated using the procedure detailed in section 2.3 with polynomial orders N = 3 and N = 6, respectively. All had the same initial position t in = 0 and angle θ in = 0. The colour coding represents each trajectory's likelihood. The right panel in figure 4 is essentially a cross-section through the left plot at depth u = 10 cm where each green dot corresponds to one trajectory and its likelihood is represented along the y-axis. The red histogram was obtained by binning the transverse positions of the trajectories. Note that for better visibility, only 10 4 trajectories are shown as dots while position information from 10 5 trajectories were used to generate the histograms. The blue histogram was obtained based on Monte Carlo simulated proton trajectories such as the ones shown in figure 3 left. The widths of the dashed black Gaussian distributions are calculated analytically according to equation 25.

Figure 5 shows polynomial coefficients obtained from fits to Monte Carlo simulated proton trajectories (see section 2.8) and from analytical sampling (see section 2.3). A polynomial trajectory was fitted to each simulated trajectory as shown in figure 3 yielding a coefficient vector a. Each dot in the scatter plots corresponds to a pair of components of one vector a, i.e., of one simulated trajectory. The histograms in the upper row show marginalised distributions of one coefficient only. It is interesting to observe that these latter are appreciably broader for the Monte Carlo data than for the analytically generated polynomials, i.e. a 2 n is larger, while the spatial distributions shown in figure 4 match well. The reason is that the coefficients from Monte Carlo are also more anti-correlated ( a n a n < 0) than the analytical ones which lowers the spatial variation t 2 (u) (see equation 10).

RESULTS

Figure 4: Left: Analytically generated proton trajectories (section 2.3) using a polynomial order N = 3 (upper row) and N = 6 (lower row). The entry energy was assumed to be 200 MeV. Right: distribution of transverse positions at 10 cm depth obtained from analytical trajectories and those simulated with Monte Carlo under identical conditions. Each green dot corresponds to one trajectory in the left panel.

RESULTS

Figure 5: Pair-wise distribution of polynomial coefficients and marginalised distributions (upper row). Blue dots were obtained by fitting polynomials of order N = 6 to Monte Carlo simulated trajectories (see figure 3). Red dots were generated analytically (see section 2.3) using the same polynomial degree.

RESULTS

Figure 6:

RMSE between simulated proton trajectories and different MLP estimates as a function of depth in 20 cm water phantom for 180 MeV (left) and 200 MeV (right) protons. For the yellow, grey, and red curves, we used our polynomial expression (equation 21) and for the blue dash-dotted one, we used the expression from [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] integrating the scattering power T (u) numerically without any polynomial approximation.

MLP estimation in proton CT reconstruction

Figure 6 shows the RMSE between proton trajectories simulated with Monte Carlo and several MLP estimates: on one hand using a polynomial approximation to the scattering power T (u) (equation 17) and consequently the polynomial expression for the MLP (equation 18) for different orders M ; and on the other hand using the general MLP expression (equation 14) and integrating the T (u) function numerically. The boundary conditions t in , t out , θ in , and θ out , were obtained from the Monte Carlo data. A total of 10 5 trajectories were simulated and for each of them the MLP was estimated. For M = 5 (not shown), the RMSE is indistinguishable from the one obtained when integrating T (u) numerically. Likewise, the RMSE obtained with the conventional MLP (not shown) is identical to the results for M = 5. For 200 MeV protons (right panel), already M = 1 yields a RMSE as good as for higher orders. For 180 MeV (left panel) on the other hand, the RMSE remains slightly higher even for M = 1 or M = 2 because the polynomial approximation of T (u) is less accurate in this case (see figure 2).

Reconstructed proton CT images

Figure 7 shows reconstruction results of the Catphan phantom obtained with different MLP estimates at 200 MeV and 180 MeV. On one hand, we used the expression given in equation 21 and approximated the scattering power T (u) function with polynomials of varying degrees (see figure 2). On the other hand, we estimated the MLP through the expression from [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] approximating T (u) by a fifth order polynomial. The root mean squared difference between the proton CT images using our polynomial MLP estimate and Schulte's is about 3.2 % for N = 5 and 4.5 % for N = 0, compared to the noise level of the images of about 3.3 % (upper row). A faint ring-shaped structure is visible in the difference image between N = 5 and N = 0 (upper right) where the inserts are located, which is to be attributed to the slightly more blurry line patterns in the latter case. This is coherent with the detail images in the centre row which are a bit more blurred for N = 0 than when using higher order approximations of T (u). At 200 MeV (lower row), the line patterns appear almost identical in all cases.

Reconstruction was performed on a MacBook Pro (2017) on a single Intel i7 processor. MLP estimation took 2.7 s per projection on average (about 68000 protons per projection) using Schulte's expression. This amounts to 32 min for the 720 projections. Evaluating our polynomial expression with order 1 took 0.1 s per projection on average, or 72 s for all 720 projections. This 4. DISCUSSION is essentially negligible compared to the time needed for binning data into projection images, reconstructing, and reading/writing data. The gain in computation speed is more than what would be expected by a pure comparison of floating point operations probably because the compiler is able to optimize the code more efficiently given the simpler structure of the mathematical operations (see section 2.7). We underline that similar optimisation may be possible for Schulte's expression and the reported computation speeds should serve as indications only.

DISCUSSION

The likelihood functional introduced in this work provides a way to generate proton trajectories in medium as polynomial functions as explained in section 2.3. The medium is required to be homogeneous in the direction transverse to the beam, but may show heterogeneities along the beam, as we will discuss further down. Each polynomial proton trajectory is represented by an N + 1 dimensional coefficient vector, where N is the polynomial order. The first two coefficients represent the entry positions and angle of the proton. The method is computationally very efficient because it only requires sampling the coefficient vectors according to a normal distribution followed by a multiplication with an eigenvector matrix for each coefficient vector. A limitation of the method is that nuclear interactions are neglected. In this sense, analytically generated projections correspond to Geant4 Monte Carlo, e.g., with the emstandard physics list, i.e., electromagnetic interactions only.

One potential application of the trajectory sampling method could be the analytical generation of proton CT data without running a full Monte Carlo simulation. Projection images (containing water equivalent thickness) could be obtained by integrating relative stopping power maps along trajectories by ray tracing methods. The technical effort necessary to implement such an analytical proton CT simulation as in [START_REF] Gianoli | Analytical simulator of proton radiography and tomography for different detector configurations[END_REF] was beyond the scope of this work. Another application of the analytical trajectory generation could be in Monte Carlo treatment planning systems. The most abundant process when transporting a proton through a patient geometry is MCS. Instead of performing this explicitly by successively deflecting the proton path by small angles, the overall MCS contribution to the trajectory could be sampled analytically. Discrete nuclear events would need to be simulated additionally and protons towards their end of range would potentially need to be transported explicitly. Developing such methods was beyond the scope of this work, but would deserve further investigation in our opinion.

Figure 4 confirms that the spatial distribution of analytically sampled trajectories resembles closely the one generated by Monte Carlo. Interestingly, if the polynomial degree N is low, such as in the upper right panel where N = 3, the generated trajectories are intrinsically constrained to be relatively smooth and tend to have a high likelihood especially if concentrated around the central axis of the beam. On the other hand, with N = 6 (lower right panel), the polynomials have more freedom to generate small scale curvature so that their likelihood is lower on average and less correlated with the position relative to the beam centre. Furthermore, figure 5 shows that the distribution of coefficients obtained by fitting polynomials to Monte Carlo trajectories resembles closely the analytical (N -1)-dimensional Gaussian distributions. The higher order coefficients (a 5 , a 6 ) obtained from Monte Carlo spread to slightly larger values than the analytically generated ones. We note in this context that the angular distribution due to MCS is actually not purely Gaussian, but contains a Rutherford-like component leading to higher tails [START_REF] Scott | The Theory of Small-Angle Multiple Scattering of Fast Charged Particles[END_REF] which are expected to generate more small scale curvature (larger high order coefficients). It is also known that the MCS models implemented in Geant4 (and other Monte Carlo codes) bear some uncertainty compared to theoretical MCS models [START_REF] Makarova | Comparison of Geant4 multiple Coulomb scattering models with theory for radiotherapy protons[END_REF]. In any case, the parametrisation of the scattering power T (u) could be refined to optimise the match between analytical and Monte Carlo trajectories. For this work, we aligned the definition of T (u) with the ones used in other contributions on MLP estimation.

We derived an expression for the MLP by maximising a likelihood functional using the Euler-

CONCLUSION

Lagrange formalism. The resulting general expression (equation 14) is strictly identical to the conventional ones in [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF] and [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] (equation 23), except for the logarithmic pre-factor in Schulte's expression for the integrated scattering power (see section A).

In fact, it would be possible to derive our MLP starting from the conventional equation 23.

When the function T (u) is approximated by a polynomial, which is common practice in proton CT reconstruction, the MLP expression automatically becomes a polynomial (equation 18). The order M of the T (u) polynomial determines the order N = M + 3 of the MLP polynomial.

For M = 5, as used by [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF] and [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF], the proton CT images of the Catphan phantom reconstructed with our MLP and with the conventional expression are virtually identical (figure 7). The residual noise is mainly due to the additional pre-factors in Schulte's MLP formula which are not present in ours (see equations 25 to 27). There are also slight numerical differences in the polynomial fits to the 1/β 2 p 2 (u) function used in Schulte's formalism and to the scattering power T (u) used in our method, respectively. A cubic MLP, such as used by [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF], is equivalent to M = 0 and thus implies neglecting energy loss in the medium. Inside the medium, the cubic MLP expression yields a less accurate estimate than an MLP estimate with M > 0 (figure 6). Accordingly, the reconstructed Catphan images have a lower spatial resolution (visually assessed), although only slightly. This is coherent with the results reported by [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF]. [START_REF] Collins-Fekete | Developing a phenomenological model of the proton trajectory within a heterogeneous medium required for proton imaging[END_REF] also estimate the MLP with a cubic polynomial, however only impose the constraints on entry and exit position, t MLP (u = 0) = t in and t MLP (u = u out ) = t out . The propagation angles at θ MLP (u = 0) and θ MLP (u = u out ) are scaled through an empirical Monte Carlo based calibration.

In view of the methods reported in this work, an alternative approach would be to select the desired order M and obtain N as a consequence rather than fixing it a priori. The choice of M generally depends on the water equivalent thickness of the object compared to the proton range at the beam energy used for image acquisition (see figure 2). Our results (see figure 7) suggest that M = 0 would be sufficient for a 20 cm (water equivalent) thick object if the proton range is 26 cm (200 MeV), but that M = 1 (or higher) would be preferable when the range is only 21.7 cm (180 MeV). We note that ideal trackers were used in the simulation and that the measurement uncertainty of a real proton CT scanner would additionally degrade the image resolution [START_REF] Krah | A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution[END_REF].

Furthermore, we evaluated the scattering power T (u) only in water in this study so that the analytically generated trajectories as well as the MLP formulae are strictly valid only in water. The same is true in conventional MLP estimation methods used in proton CT reconstruction. We underline, however, that material heterogeneities along the beam path, i.e. which are mainly a function depth u, could be easily considered. An example could be protons passing trough different slabs of soft tissue and bone. To this end, the protons' energy loss would need to be calculated in the heterogeneous geometry to adjust the β 2 p 2 -term in T (u) and the X 0 (u) factor would need to reflect the varying radiation length of the material, similar to what has been suggested by [START_REF] Collins-Fekete | Extension of the FermiâĂŞEyges most-likely path in heterogeneous medium with prior knowledge information[END_REF]. Such an approach would provide a way to sample proton trajectories and estimate an MLP in media with the above described heterogeneities. Clearly, this would require some prior knowledge about the imaged object.

If variation of material properties in depth are considered by appropriately adjusting the β 2 p 2 -term in T (u) and the radiation length X 0 (u), the polynomial degree M will likely depend on the degree of heterogeneity. The impact of material heterogeneities on the MLP has been studied elsewhere [START_REF] Wong | The Effect of Tissue Inhomogeneities on the Accuracy of Proton Path Reconstruction for Proton Computed Tomography[END_REF][START_REF] Khellaf | Effects of transverse heterogeneities on the most likely path of protons[END_REF] and was not the subject of this work.

CONCLUSION

We proposed a functional to quantify the likelihood of proton trajectories under the influence of multiple Coulomb scattering and achieved two things: first, the likelihood functional provides an efficient way to analytically sample proton trajectories in homogeneous media instead of gen-erating them trough a Monte Carlo simulation. Second, maximising the functional yields an expression for the most likely path of a proton which can be used in proton computed tomography reconstruction. When the protons' energy loss is approximated by a polynomial function of depth, the most likely path is a polynomial itself. Our most likely path expression is analytically equivalent to the conventional ones available in the literature. Our results suggest that it is sufficient to describe the energy loss with a linear or at most quadratic function of depth instead of a fifth order polynomial. Treating the proton energy as constant leads to a slightly degraded spatial resolution. Evaluating our most likely path expression requires fewer floating point operations than the conventional ones.

Figure 1 :

 1 Figure 1: Depth dependent inverse scattering power 1/T (u) = β 2 (u)p 2 (u)X 0 (u)/E 2 0 for 200 MeV protons in a homogeneous water phantom. Coloured curves are polynomial fits.

Figure 2 :

 2 Figure 2: Depth dependent scattering power T (u) = E 2 0 /β 2 (u)p 2 (u)X 0 (u) for 200 MeV (upper) and 180 MeV (lower) protons. The dashed vertical line on the right side indicates the proton range at the respective energy.

Figure 7 :

 7 Figure 7: Upper row: Simulated proton CT image (left) of the Catphan phantom at 180 MeV energy reconstructed using equation 21 and a fifth order polynomial to fit the scattering power T (u); difference to the proton CT image obtained with the MLP expression from Schulte et al. (2008) (centre) and using an order M = 0 polynomial to approximate T (u) (right). Center (180 MeV) and lower (200 MeV) row: Details of the phantom indicated by red and yellow boxes in the upper left panel; intensity profiles across the two line patterns averaged along the line direction.

Acknowledgements

The work of Nils Krah was supported by funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 753370. This work was performed within the framework of the project DIC20161236452 of the "Fondation pour la Recherche Médicale" (FRM), of the SIRIC LYriCAN Grant INCa_INSERM_DGOS_12563 and of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the ANR.

A Relation with conventional MLP expression

We summarize the steps to link our MLP expression (equation 14) with the ones reported by [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF] and [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF], which we state here for completeness:

where

The matrices R 0 and R 1 are essentially small angle rotation matrices and Σ 1 and Σ 2 quantify the spatial and angular dispersion in a certain depth due to MCS. Their components are calculated through the following integral expressions:

One obtains the components of Σ 2 by replacing u in with u and u with u out in the equations above. We recall that we have chosen u in = 0 for simplicity in our work. The quadratic prefactor containing the logarithmic term is reported only by [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF], not by [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF].

To derive the relation between our MLP formulae and the conventional one, we observe that the integral expressions in equations 13 and 14 can be linked to the conventional ones (equations 25, 26, 27; neglecting the logarithmic pre-factor) as