P Schmidt-Wellenburg
email: philipp.schmidt-wellenburg@psi.ch

J Bossy

E Farhi

M Fertl

K K H Leung

‡ A Rahli

T Soldner
email: soldner@ill.eu

O Zimmer

Experimental study of ultracold neutron production in pressurized superfluid helium

Keywords: Ultracold neutron production, pressurized superfluid helium, source of ultracold neutrons

We have investigated experimentally the pressure dependence of the production of ultracold neutrons (UCN) in superfluid helium in the range from saturated vapor pressure to 20 bar. A neutron velocity selector allowed the separation of underlying single-phonon and multiphonon processes by varying the incident cold neutron (CN) wavelength in the range from 3.5 to 10 Å. The predicted pressure dependence of UCN production derived from inelastic neutron scattering data was confirmed for the single-phonon excitation. For multiphonon based UCN production we found no significant dependence on pressure whereas calculations from inelastic neutron scattering data predict an increase of 43(6) % at 20 bar relative to saturated vapor pressure. From our data we conclude that applying pressure to superfluid helium does not increase the overall UCN production rate at a typical CN guide.

I. INTRODUCTION

Ultracold neutrons (UCN) have energies below 300 neV and can be stored for long observation time in magnetic or material bottles [START_REF] Golub | Ultra-Cold Neutrons[END_REF]. This peculiarity makes them attractive for high precision measurements of fundamental properties of the neutron which are relevant for particle physics and cosmology [START_REF] Dubbers | [END_REF]. The most prominent example is their use for searches of a neutron electric dipole moment [3][4][5][6][7][8] which would give direct evidence of CP-violation beyond the Standard Model of particle physics. Other projects aim to improve the knowledge of the neutron β-decay lifetime which is crucial for calculations of big bang nucleosynthesis [9]. Together with correlation measurements in neutron decay it contributes to a more precise understanding of weak semileptonic processes in the first quark generation [10,11].

New developments might follow from observations and manipulations of quantum states of UCN in the Earth's gravitational field [12][13][14]. All of these endeavors greatly benefit from an increase of UCN counting statistics in their experiments.

Current projects to increase the density of UCN for physics experiments employ neutron converters of superfluid helium (He-II) [15][16][17][18][19][20] and solid deuterium [21][22][23][24][25][26]. They exploit the concept of a superthermal source [27] where the converter temperature can be much higher than the UCN "temperature". Upscattering of UCN is suppressed by the Boltzmann factor if the energy of excitations in the converter is large compared to its temperature. Using liquid He as converter was proposed in [28] and first experimentally realized in [29].

In He-II the main contribution to UCN production is due to single-phonon excitations (see [17] for a separation of UCN production in single-and multiphonon contributions), which occur at the crossing point λ * of the dispersion relations of the free neutron and of He-II.

At saturated vapor pressure (SVP), λ * = 8.92(2) Å (derived from the data in [30], [50]).

Application of pressure shifts λ * to lower values, where the differential flux dφ/dλ is higher for a typical neutron beam in a guide coupled to a liquid deuterium CN source, which typically has its flux maximum around 4 Å. Combined with an increase in density of He-II with pressure one may anticipate an overall gain in UCN production.

However, a quantitative analysis based on a method described in Ref. [31] using inelastic neutron scattering data at SVP and 20 bar from [30] predicted a net decrease in single-phonon UCN production [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF]. On the other hand, these same calculations also predicted an increase in the production from multiphonon excitations at 20 bar. We wished to test these calculations and quantify the UCN production in pressurized He-II by a direct measurement. Independent of UCN production, another motivation for performing these studies is the larger dielectric strength of pressurized He-II [START_REF] Hara | [END_REF]. This is of particular interest for searches of the neutron electric dipole moment within He-II using electrical fields of several 10 kV/cm [35].

II. EXPERIMENTAL SETUP

The experiment took place at the CN beam facility PF1b [36] of the Institut Laue-Langevin, France. A schematic diagram of the experimental setup is shown in Fig. 1. The setup for the He-II UCN production target and UCN extraction was similar to the one described in [18][19][20]: The target was located inside a cryostat equipped with a commercial two-stage Gifford McMahon cold head with a cooling power of 1.5 W at 4.2 K. In the cryostat, gaseous He with 99.999% purity was liquefied and cooled below the λ-transition to superfluidity by a continuous 4 He evaporation stage. 3 He was removed by passing the superfluid through a superleak (for technical details on the used design, see [19]) into the UCN production volume. Ultracold neutrons were extracted by an inverted U-shaped stainless steel UCN guide to a 3 He UCN detector. For further details, see [18,20,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF].

In order to allow for measurements at high pressure, a dedicated He-II container was used as converter vessel. It consisted of a stainless steel tube (inner diameter 66 mm, length 193.8 mm) with two spherical end caps from aluminum (radius ∼ 33 mm, wall thickness 3 mm) for passage of the CN beam. Its total volume was 810 cm 3 and the intersection volume with the neutron beam about 180 cm 3 . The cylinder was coated with natural nickel and the end caps had nickel half-spherical shells inserted on the inside to increase the neutron optical potential for UCN storage. The UCN flapper valve used in [18] for window-less extraction was replaced by an aluminum window (0.1 mm thickness) supported by a stainless steel disk with an array of holes (effective hole area 127 mm 3 ). The 3 He evaporation stage serving to further cool the production volume had to be operated with 4 He due to the higher heat load and a beam monitor with circular orifices d) before it entered the cryostat through a series of thin aluminum windows e). The He-II target was kept in the Ni coated stainless steel pressure container g); CN entrance and exit sides from aluminum were covered with nickel half-spheres q). Ultracold neutrons were guided by a polished stainless steel tube f), which was separated by an aluminum window from the production volume, to the 3 He UCN detector m). The CN beam exited the cryostat through aluminum windows h) before entering a chopper k) and TOF-detector l) for time-of-flight spectral analysis. Pressure was applied via the filling line p) connected on the bottom of the pressure container and measured on a static line with a gauge o).

of the high pressure volume. Therefore the target temperature reached only 1.1 K instead of 0.7 K obtained in [18]. The temperature of the He-II was measured using a calibrated Cernox [START_REF]Model: CX-1030-CU[END_REF] sensor screwed tightly to the outside of the He-II vessel. The pressure was applied from the gaseous He supply via the liquefaction chain through the superleak and measured with a room-temperature pressure gauge calibrated for p < 100 bar relative to atmospheric pressure. The gauge was connected to the converter via a thin stainless steel tube which could also be used for pumping to establish SVP conditions.

A neutron velocity selector [37] was used to scan the incident CN wavelength range. It had to be used in two orientations, -5 • and +5 • with respect to the incident beam, as the normal position (0 • ) did not give access to the entire wavelength range of interest. Hence, two scans were performed: from 3.5 Å to 6.5 Å and from 5.5 Å to 10 Å. Depending on the angular and rotational settings of the velocity selector a wavelength spread of ∆λ(FWHM) λ ≈ 0.14, 0.08 was obtained for -5 of the converter (τ = 2.92 (7) s for SVP at T = 1.08(2) K). The beam-independent background r 0 = 0 +3.9 -0.0 × 10 -3 s -1 was determined in dedicated measurements with CN-shutter closed. This parameter was kept fixed in the fits. It was found that the storage lifetime was independent of pressure (see Fig. 3b). Increasing the temperature to 1.48 K decreased τ to 1.7(1) s. This temperature behavior combined with results from UCN Monte Carlo simulations with Geant4UCN [38] and simple gas kinetic arguments showed that the storage lifetime was dominated by UCN escaping through the extraction holes. Hence, changes in count rate due to the small temperature fluctuations of ±0.02 K were negligible.

In phase (ii), the background was increased by a beam-dependent component r beam causing a constant offset to the UCN build-up prompt with opening the CN shutter. A systematic error will follow if this offset is not properly taken into account. The component r beam was determined by fitting the build-up with an exponential with constant background r = r 0 + r beam . It was expected to be proportional to the CN monitor count rate since CN monitor and UCN detector both had a low efficiency proportional to 1/v for CN. Therefore r beam was determined for four different velocity selector settings (maximum of single-phonon production and three settings with low ratio of UCN production to CN monitor count rate).

A fit of the four values with r beam ( ṅCN ) = k ṅCN , where ṅCN was the CN count rate of the beam monitor (12000 s -1 < ṅCN < 17000 s -1 ), yielded k = (5.6 +6.0 -5.6 ) × 10 -7 . All data of phase (ii) was corrected for both background contributions, by typically 0.01 s -1 .

In addition to these scans the UCN production rates across the single-phonon peak at SVP were measured several times during the entire experiment to check the reproducibility.

The UCN count rate was found to decrease with time. This is suspected to be due to gas freezing out on the thin aluminum window at the extraction holes or on the surface of the extraction guide. An exponential decrease with respect to the start of the first measurement (t = 0): N (t) = N 0 exp -t T + A 0 with a large offset A 0 = 1.95(10) × N 0 and time constant T = 104 (12) h was found to describe this effect well and thus was used to correct the data.

The corrected UCN counts were normalized to the particle flux as measured by the beam monitor (monitor count rate is corrected for dead time and 1/v efficiency). Time-of-flight measurements (see Fig. 2) were used to determine a central wavelength for each setting of the velocity selector. The measured TOF spectra were deconvoluted taking into account the opening function of the chopper and the finite thickness of the TOF detector. A logarithmic normal distribution was found to describe the skewed deconvoluted spectra better than a symmetric Gaussian and was used to find the position of the maximum of each spectrum.

Hence, each UCN production rate, which is in fact an integral UCN production rate of the specific CN spectrum of this velocity selector setting, was assigned to one wavelength λ. The uncertainty of the wavelength measurement was estimated as ∆λ λ = ±1 %.

IV. RESULTS

The measured wavelength-dependent UCN rates for all pressures, normalized to the in- gets less pronounced with pressure at the available wavelength resolution. Note that the wavelength resolution was limited by the resolution of the velocity selector (the same for the range of all single-phonon peaks shown). With pressure the single-phonon peak decreases in intensity and moves towards shorter wavelengths.

cident
The area under the peak is proportional to the UCN production rate by single-phonon excitation. The apparent peak width comes from the convolution of the incident CN spectra with the extremely narrow single-phonon scattering function (the linewidth of the roton excitation in 4 He at ∼1 K is of the order of 1 µeV, see e.g. Ref. [39], compared to the resolution of ∼ 230 µeV of the velocity selector at 8 Å). The peak was fitted with a logarithmic normal distribution in q-space. This ansatz fitted best the CN spectra and was used here again since the single-phonon UCN production rate is proportional to the CN intensity at λ * for each incident CN spectrum. Note that fits with logarithmic normal distributions gave for all pressures reduced χ 2 values closer to unity than a symmetric Gaussian distribution. The fitting range was restricted to wavelengths larger than 7.5 Å, hence excluding contributions from multiphonon processes. The positions of the maximum of each fit function were observed for all pressures at values smaller than expected from the crossing point of the dispersion relations of He-II and the free neutron. They are listed in Tab. I together with the UCN production rates normalized to the rate at SVP. All single phonon UCN production rates R I were taken from the integrals over the fitted curves for λ > 7.5 Å.

The obtained fits were then subtracted from the data in order to extract the UCN production rates due to multiphonon processes. These are shown scaled up in Fig. 5 and found to exhibit a broad, pressure-independent maximum at approximately 6.25 Å. For estimating the relative change of UCN production due to multiphonon excitation and multiple scattering with pressure we used a simple sum over data points for λ < 7.5 Å, which was found to remain constant within errors in this range (see Table I).

V. DISCUSSION

In our experiment we observed the anticipated decrease of UCN production R I due to single-phonon processes with pressure. This behavior can be well understood by taking a look at

R I ∝ N βλ * 4 S * dφ dλ λ * , ( 1
)
as derived in an unpublished note of Pendlebury (see also [START_REF] Golub | Ultra-Cold Neutrons[END_REF]), where the single-phonon scattering intensity S * = S(q * , ω * ) decreases from S * (SVP) = 0.118 (8) to S * (20 bar) = 0.066( 6) [32]. This change of the structure factor explains most of the observed decrease.

Smaller contributions come from the change of λ * and from reduced overlap between the dispersion relations of He-II and of the free neutron taken into account by the parameter β(SVP) = 1.42 to β(20 bar) = 1.21 (values deduced from the scattering data [30], see [START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF]).

The increase in helium number density N [10 22 cm -3 ] from N (SVP) = 2.1836 to N (20 bar) = 2.5349 (calculated from the mass densities at 0.1 K from [40]) cannot compensate these effects. Note that the factor dφ dλ λ * does not contribute to the normalized data in Fig. 4. Its inclusion for a typical CN spectrum does also not compensate the discussed decreases. This behavior confirms the results from calculations, which agree within errors with the experimental results, see Tab. I and Fig. 6.

We observe a displacement of all measured single-phonon production rate peaks to lower wavelengths than anticipated from the crossing points of the free neutron dipsersion curve with the one of He-II. This can be explained by an unintentional displacement of TOF monitor by 5 cm, the width of the detector. This most probably has happened during the setup for the second measurement period (the only one which delivered data). For this reason we present all our results in table I using a "corrected" wavelength scale λ cor = 1.025•λ.

In the multiphonon range we observed a discrepancy between the measured and the calculated pressure dependence of the UCN production rate (see Tab. I and Fig. 5). From calculations we would have expected an increase with pressure, whereas no increase could be observed. The observed pressure independence of the storage time constant does not indicate a change of the UCN extraction probability from the converter vessel with pressure. Note that our experimental approach provides a direct relative comparison of the production rates at different pressures: For a fixed wavelength, the same incident beam was used. Pressuredependent UCN losses (e.g. due to cracks in the coating of the storage vessel) are independent experiment from Ref. [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF] (see text), and UCN production rates compared with results given in Ref. [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF] of the CN wavelength.

p λ * λ * cor s-ph m-ph p λ * s-ph m-ph [bar] [Å] [Å] λ < 7.5 Å [bar] [Å] λ < 7.
One explanation for the discrepancy in the multiphonon region could be a change in the dynamic structure factor S(q, E) with temperature. The calculations Ref. [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF] employ scattering data measured at 0.5 K whereas our measurements were done at 1.1 K. Increasing temperature results, for example, in a broadening of the single-phonon excitation. This increase is small at low temperatures but becomes large towards the λ-transition [30]. In [30] the temperature dependence of the single-phonon excitation linewidth for different pressures was determined from fits to the measured dynamic structure factor. No significant change of these values was found below about 1.3 K. Likewise, the multiphonon scattering varies only weakly with temperature [30]. Therefore we consider the scattering data measured at 0.5 K as adequate to predict UCN production at 1.1 K. Note that temperature also influences UCN from this experiment ( ) and from calculation from scattering data (×) [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF]. Both dependencies have been normalized to the respective production rate at SVP. up-scattering but this effect cannot explain the observed discrepancy as it is independent on the CN wavelength, and would also shift the single-phonon data from expectation.

Assuming no significant change of the dynamic structure factor between 0.5 K and 1.1 K, the observed discrepancy can only be explained by an effect depending on both, the CN wavelength and the He pressure, such as inelastic scattering to other than UCN energy.

Inelastic scattering has two types of consequences: (i) reduction of the beam intensity averaged over the production volume with respect to the incident one, as scattered neutrons may leave the production volume sideways; and (ii) change in CN energy which may influence the probability of UCN production for the scattered neutron (multiple scattering). Both effects are not taken into account in the calculations.

In order to estimate the reduction of the beam intensity (i) we use the wavelength-and temperature-dependent transmission of neutrons through liquid He at SVP measured by Sommers et al. [41]. We calculate the volume-average CN flux in the multiphonon range (from 3.5 Å to 7.5 Å) inside the production volume relative to the incoming flux:

t m-ph = m-ph dφ dλ t(λ) dλ m-ph dφ dλ dλ . ( 2
)
Values between 3.5 Å and 4.52 Å were obtained from linear extrapolation of the data given there. We obtain t m-ph (SVP) = 0.89. For 8.9 Å, the correction is t s-ph (SVP) = 0.98. Note that for a longer production volume, the correction would be larger. We define multiphonon and single-phonon production rates, R II and R I , as integrals over the respective regions in the CN-normalized UCN production spectrum shown in Fig. 4. In order to compare the observed ratio R II /R I SVP = 0.14 [START_REF] Golub | Ultra-Cold Neutrons[END_REF] with the calculations, we need to correct the prediction with the geometry specific attenuation: [START_REF] Golub | Ultra-Cold Neutrons[END_REF].

t m-ph /t s-ph • R calc II /R calc I SVP = 0.11
A reduction of CN transmission t m-ph with increasing pressure could explain the observed discrepancy between the measured and the calculated pressure dependence of multiphonon UCN production. Unfortunately, a quantitative confirmation of this hypothesis is not possible, due to lacking transmission data for higher pressures. Also scattering data for different pressures [30] do not cover the high-q range needed to calculate the transmission via total scattering.

The change of the CN energy by inelastic scattering (ii) depends on the incident wavelength and the pressure. Scattering to the wavelength λ * is kinematically possible and would increase the probability of UCN production by the scattered CN; however, the width of singlephonon UCN production is very small. Multiple scattering was not taken into account in the calculations [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF]. Its contribution to the UCN production in our setup was estimated by Monte Carlo simulations [42] based on the McStas package [43,44]. For SVP, multiple scattered CN were found to contribute about 6 % to the UCN production rate in the multiphonon region. This fraction decreases slightly with increasing pressure, to 5 % for p ≥ 10 bar. The weak pressure dependence of this small contribution cannot explain the discrepancy between observed and expected pressure dependence of multiphonon UCN production. It changes the calculated relative production rates at SVP t m-ph /t s-ph • R calc II /R calc I SVP to 0.12(1) .

However, it should be noted that the accuracy of the simulations and calculations is limited by the available scattering data: Single data sets cover only a limited (q, ω)-range or provide a limited resolution. Therefore, data sets [START_REF] Bossy | Data from experiment TEST-1715 at IN5 (incident wavelength 3 Å) at the ILL[END_REF][START_REF] Bossy | Data from experiment 6-01-206 at IN5 (incident wavelength 5 Å) at the ILL[END_REF][START_REF] Gibbs | Data from experiment 6-01-140 at IN6 (incident wavelength 4.1 Å) at the ILL[END_REF] from different instruments and with different resolution settings had to be combined in the simulations [42]. The scattering data [30] used for the calculations [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF] were obtained with two different incident wavelengths at IN6: 4.1 Å for SVP and 4.6 Å for higher pressures. Possibly uncorrected effects may contribute to the observed discrepancies.

VI. CONCLUSION

We have investigated CN-wavelength-dependent UCN production in pressurized He-II for pressures up to 20 bar. The single-phonon UCN production peak moves to shorter wavelength and decreases in intensity, whereas the absolute contribution from multiphonon processes stays constant within the experimental uncertainties. We have found that predictions from calculations based on inelastic neutron scattering data [32,[START_REF] Schmidt-Wellenburg | Production of Ultracold Neutrons in Superfluid Helium Under Pressure[END_REF] agree well with the behavior of the single-phonon UCN production rate. However, the same calculations disagree with our measurement in the multiphonon range. There is no indication for a significant change in the dynamic structure factor between 0.5 K and 1.1 K that could invalidate the calculations for our experiment. However, the disagreement could be caused by a pressure-and wavelength-dependent attenuation of the CN beam in the He-II. This hypothesis could not be tested due to a lack of transmission data at higher He pressure. Note that the ratio of multiphonon to single-phonon UCN production would reduce for longer production volumes.

The calculations also did not take into account modifications of the CN spectrum by inelastic scattering inside the He-II, but Monte-Carlo simulations showed that the contribution of UCN production by multiple scattered CN and its pressure dependence are too small to explain the observed disagreement.

The disagreement of the position of the single-phonon peak at SVP with earlier measurements by Yoshiki and coworkers [START_REF] Yoshiki | [END_REF] and Baker and coworkers [17] can most probably be explained by an accidental displacement of the TOF detector in our experiment.

In general, the relative UCN production rate due to single-phonon and multiphonon processes depends upon the incident CN spectrum. For the spectrum in our experiment we have found that (30 ± 2) % of the total UCN production over all incident wavelengths are from multiphonon processes. Baker and coworkers [17] found a multiphonon contribution of (24 ± 2)% for a different CN beam. Note that, according to our previous discussion, also the length of the UCN production volume (which was 326 mm in the experiment of Baker and coworkers) influences this ratio.

Using the entire spectrum of a white beam might still be advantageous for a powerful UCN source, as monochromatization by a crystal not only cuts out a narrow wavelength band around 8.92 Å, but also reduces the intensity of the Bragg reflected beam substantially (factor 2-5 [20,32]). In whatever case, monochromatic or white incident beam, best results are obtained for He-II at SVP; without even taking into account the increase in technical complexity when going to higher pressures. For future searches of the neutron electric dipole moment it might still be beneficial to apply a pressure of a few 100 mbar to increase the dielectric strength without loosing significantly in UCN production. A future perspective to increase UCN production might be the use of solid helium as UCN converter, which was not possible to investigate during this experimental run.

Figure 1 :

 1 Figure 1: (Color online) Schematic of the setup (not to scale): The CN beam passed through the velocity selector a), a secondary CN guide b) which could be closed off with a fast shutter c),

  beam monitor was used to open and close the beam (opening and closing time less than 1 s). Before entering the beam monitor the beam's cross section was reduced to ∅ = 30 mm with two circular apertures made from B 4 C. Images of the beam at various distances behind the monitor were taken before the cryostat was set up. They showed that the beam passed through the production volume without touching its cylindrical sidewall for all selected incident wavelengths. A time-of-flight (TOF) system for spectral analysis of the incident spectrum was set up behind the cryostat. It consisted of a chopper -a disc rotating at 1000 RPM with four radial rectangular slits (width 3 mm) passing by a fixed aperture of the same dimensions -and a 3 He CN-monitor with a detection efficiency proportional to 1/v (v is the neutron velocity) placed 2 m downstream. The corresponding wavelength spectra for some selected velocity selector settings are displayed in Fig.2.III. MEASUREMENTSScans of UCN count rate versusCN wavelength were performed at absolute pressures and temperatures (p [bar], T [K]) = (6, 1.105), (11, 1.106), (16, 1.112), (20, 1.109), and (SVP, 1.08), with temperature fluctuations ∆T ≤ 0.02 K during individual scans. The UCN production rate was first measured in the long wavelength range with the velocity selector angle set to +5 • for all pressures. Then the measurement of the UCN production rate in the short wavelength range (-5 • ) followed. For each pressure and angular setting, three updown scans were performed with the same set of velocity selector frequencies corresponding CN wavelength (Å)

Figure 2 :

 2 Figure 2: (Color online) Time-of-flight spectra of the CN beam measured with empty production volume, for various velocity selector rotation speeds at -5 • ( -) and +5 • ( ---) selector alignment.Each raw spectrum was normalized to area unity and then multiplied with the respective monitor count rate corrected for dead time and 1/v efficiency. The difference in integral counts for -5 • and +5 • was due to a larger transmission of the velocity selector at -5 • .

Figure 3 :

 3 Figure 3: (Color online) Count rate observed in the UCN detector (a). Shown is the combined data of all measurements at SVP and λ ∼ 8.7 Å (maximum of UCN production, λ cor = 8.92 Å). The vertical lines indicate the three different phases. The solid line in phase iii) is the best fit (see text) to the data with a reduced χ 2 = 1.04. The decay time constants τ for different pressure were constant within the error bars (b): The red line (-) is the weighted mean, the blue lines (-) indicate the error of the weighted mean.

  particle flux, are shown in Fig.4. Two different regions can be distinguished for all pressures: a broad distribution for short incident wavelengths from UCN production by multiphonon processes R II and a pronounced peak at long wavelengths from UCN production from single-phonon excitation R I in He-II. These regions are separated by a minimum which

Figure 4 :

 4 Figure 4: (Color online) Measured wavelength resolved UCN production rate per incident CN for p = SVP( ), 6 ( ), 11 ( ), 16 ( ), 20 bar ( ). The lines are best fits of logarithmic normaldistributions to the data with a reduced χ 2 = 1.3, 1.2, 2.0, 1.9, 2.6, respectively.

Figure 5 :

 5 Figure 5: (Color online) Ultracold neutron production per incident CN from multiphonon processes for p = SVP( ), 6 ( ), 11 ( ), 16 ( ), 20 bar ( ). Measured data (a) after subtracting the singlephonon part (see text). Calculated (b) from inelastic scattering data for p = SVP( ), 5 ( ), 10 ( ), 15 ( ), and 20 bar ( ) .

  calculated from scattering data. The UCN production rates are normalized to the incident cold neutron flux. The single-phonon (s-ph) and multiphonon (m-ph) columns show the production rates at pressure p normalized to the respective SVP results. The quoted errors of λ * from this experiment only include the uncertainties from the fits with the logarithmic normal distribution. The quoted errors of the measured UCN production rates are statistical uncertainties; in the single-phonon case uncertainties of the fit were increased to yield a reduced χ 2 of 1. The errors of the calculated values are the quadratic sums of all contributions.

Figure 6 :

 6 Figure 6: (Color online) Comparison of the pressure dependence of single-phonon UCN production

Table I :

 I Pressure dependence of measured single-phonon peak position λ * , corrected values λ *

				5 Å
	SVP 8.70(1) 8.92(set) 1	1	SVP 8.92(2) 1	1
	6 8.31(1) 8.52(1) 0.65(4) 0.93(4)	5 8.54(4) 0.75(11) 1.09(3)
	11 8.18(1) 8.39(1) 0.54(4) 0.97(4)	10 8.39(2) 0.59(10) 1.12(3)
	16 8.10(1) 8.30(1) 0.47(3) 0.97(4)	15 8.30(2) 0.47(9) 1.17(6)
	20 8.09(1) 8.29(1) 0.41(3) 1.03(4)	20 8.26(2) 0.40(8) 1.43(6)
				cor

We would like to thank the ILL staff, especially D. Berruyer, T. Brenner, P. Lachaume, P. Mutti, and P. Thomas who made this experiment possible thanks to their technical support. We would like to thank S. Mironov for his important contributions to the design of the cryogenic apparatus. This work has been funded by the German BMBF (contract number 06MT250).


Daniel Barlet

Daniel Barlet Note

Note on a differential inequality

Keywords:

come

Introduction

The motivation of this note is to look for an analogous result to the classical Gronwald lemma which gives some estimate of the growth of a positive function F near the center of a punctured disc D * of the plane which satisfies the inequality ∆F (x) ≤ C.||x|| -2 .F (x). Using the invariance by rotation of the Laplacian allows to reduce our problem to a special case of the following differential inequality

x 2 .f (x) + x.θ(x).f (x) ≤ η(x).f (x) ∀x ∈]0, 1[ (1)
for a positive function f ∈ C 2 (]0, 1[) where θ is a C 1 function in [0, 1[ and η is a continuous function on [0, 1[. We do not obtain a point-wise estimate for such an f when x goes to 0 but an integral estimate which corresponds precisely to the point-wise estimate we would wait for. But we show that such a point-wise estimate does not hold in general.

1 The estimate

The standard situation

We fix the real valued functions θ ∈ C 1 ([0, 1[) and η ∈ C 0 ([0, 1[). Let a := η(0). We shall say that we are in the standard situation when the following conditions are fulfilled 1. If q := θ(0) = 1 we assume that lim sup x→0 +

(η(x) -a).(-Log x) = α 1 < +∞.

2. If q = 1 we assume that lim sup

x→0 + (η(x) -a).(-Log x) 2 = α 2 < +∞.
We want to study the behaviour when x → 0 + of C 2 functions f :]0, 1[→ R + which are solutions of the differential inequality

x 2 .f (x) + x.θ(x).f (x) ≤ η(x).f (x) ∀x ∈]0, 1[ (1)
1.2 Reduction to the case θ (0) < 0

Assume that we are under the standard situation defined above. Then for a given real number u that we shall choose later on, define the function g :]0, 1[→ R + by letting f (x) = exp(u.x).g(x). Then f (x) = u.exp(u.x).g(x) + exp(u.x).g (x) f (x) = u 2 .exp(u.x).g(x) + 2u.exp(u.x).g (x) + exp(u.x).g (x)

so we obtain

x 2 .f (x) + θ(x).x.f (x) -η(x).f (x) = exp(u.x). x 2 .g (x) + (θ(x) + 2u.x).x.g (x) -(η(x) -u.x.θ(x) -u 2 .x 2 ).g(x)
So the differential inequality (1) for f is equivalent to the following differential inequality for the function g (which depends on the choice of u):

x 2 .g (x) + θ u (x).x.g (x) ≤ η u (x).g(x) ∀x ∈]0, 1[ (1 u ) with θ u (x) := θ(x) + 2u.x and η u (x) := η(x) -u.x.θ(x) -u 2 .x 2 .
Then we are again in the standard situation with

θ u (x) = θ(x) + 2u.x and η u (x) = η(x) + o(1/(-Log x) 2 ).
So we have θ u (0) = θ(0) and choosing u 0 we obtain that θ u (0) < 0 so the following condition will be fulfilled

• There exists ε > 0 such that θ u (x) -θ u (0) ≤ 0 on [0, ε]. (@)
Remark. The behaviour when x → 0 + of the function f and of the function g are equivalent, so to study the behaviour when x → 0 + of a solution f of the differential inequality (1) we may always assume that the condition (@) is satisfied by the function θ.

1.3 The case η(0) = 0 and q > 1 Proposition 1.3.1 We fix the functions θ ∈ C 1 ([0, 1[) and η ∈ C 0 ([0, 1[) as in the standard situation. We make the hypotheses that q := θ(0) > 1, a := η(0) = 0 and that the condition (@) is fulfilled. Then we have the following estimate:

1/2 0 t q-2 .(-Log t) -β-1 .f (t).dt < +∞ (2)
for all β > 0 satisfying β.(q -1) > α 1

Proof. First remark that the convergence in (2) is satisfied for q > 1 and any β when f is bounded near 0. This is of course the case if f is increasing on some interval ]0, c[ with c > 0. So we may assume that there exists a decreasing sequence (x ν ) ν≥0 in ]0, ε[ with limit 0 such that f (x ν ) < 0 for each ν ≥ 0 where ε is the positive number which exists from the condition (@).

Let us compute the following quantity for x ∈]0, 1/2] and β ≥ 0 fixed:

A λ (x) := 1/2 x t λ .(-Log t) -β .f (t).dt + 1/2 x t λ-1 .(-Log t) -β .θ(t).f (t).dt ≤ 1/2 x t λ-2 .(-Log t) -β .η(t).f (t).dt (1 bis)
Then, integrating by parts the first integral:

A λ (x) = 2 -λ .(Log 2) -β .f (1/2) -x λ .(-Log x) -β .f (x) + -β. 1/2 x t λ-1 .(-Log t) -β-1 .f (t).dt + 1/2 x t λ-1 .(-Log t) -β .(θ(t) -λ).f (t).dt
and integrating by parts again:

A λ (x) = 2 -λ .(Log 2) -β .f (1/2) -x λ .(-Log x) -β .f (x) -β.2 1-λ .(Log 2) -β-1 f (1/2) + β.x λ-1 .(-Log x) -β-1 .f (x) + β.(λ -1). 1/2 x t λ-2 .(-Log t) -β-1 .f (t).dt + + β.(β + 1). 1/2 x t λ-2 .(-Log t) -β-2 .f (t).dt + t λ-1 .(-Log t) -β .(θ(t) -λ).f (t) 1/2 x + -(λ -1). 1/2 x t λ-2 .(-Log t) -β .(θ(t) -λ).f (t).dt + -β. 1/2 x t λ-2 .(-Log t) -β-1 .(θ(t) -λ).f (t).dt + - 1/2 x t λ-1 .(-Log t) -β .θ (t).f (t).dt (3)
So we have:

A λ (x) = 2 -λ .(Log 2) -β .f (1/2) + (a) -x λ .(-Log x) -β .f (x) + (b) -β.2 1-λ .(Log 2) -β-1 f (1/2) + (c) β.x λ-1 .(-Log x) -β-1 .f (x) + (d) β.(λ -1). 1/2 x t λ-2 .(-Log t) -β-1 .f (t).dt + (e) β.(β + 1). 1/2 x t λ-2 .(-Log t) -β-2 .f (t).dt + (f) t λ-1 .(-Log t) -β .(θ(t) -λ).f (t) 1/2 x + (g) -(λ -1). 1/2 x t λ-2 .(-Log t) -β .(θ(t) -λ).f (t).dt + (h) -β. 1/2 x t λ-2 .(-Log t) -β-1 .(θ(t) -λ).f (t).dt + (i) - 1/2 x t λ-1 .(-Log t) -β .θ (t).f (t).dt (j)
where the term (g) is the sum of the following two terms:

2 1-λ .(Log 2) -β .(θ(1/2) -λ).f (1/2) (g 1 ) -x λ-1 .(-Log x) -β .(θ(x) -λ).f (x) (g 2 )
Now we have 1. (b) is non negative for x = x ν as f (x ν ) ≤ 0 by assumption.

2. (d) is non negative as β and f are non negative.

3. (g 2 ) is non negative for x = x ν , ν 1, when we assume that θ(x) -λ is non positive for x near enough to 0. So, under the hypothesis that the terms (b), (d) and (g 2 ) are non negative, and this is the case choosing λ = q and x := x ν < ε, the inequality (1bis) gives

1/2 xν t q-2 .(-Log t) -β-1 .f (t).Z(t).dt ≤ -2 -q .(Log 2) -β .f (1/2) + β.2 1-q .(Log 2) -β-1 f (1/2) -2 1-q (Log 2) -β .(θ(1/2) -q).f (1/2) (*)
with

Z(t) = β.(q -1) + β.(β + 1) -Log t -(θ(t) -q). (q -1).(-Log t) + β + -θ (t).t.(-Log t) -(-Log t).η(t) (**)
Now if β > sup{α 1 /(q -1), 0} we obtain that Z(t) is bounded below by a fixed positive number when t is small enough, and then (2) holds true.

Remark. Looking to the proof above, the reader may see that the following hypotheses on the function θ are enough to obtain the conclusion of the proposition, but replacing the condition β > sup{α 1 /(q -1), 0} by the (better) condition1 β > sup{α 1 /(q -1) + γ, 0}.

i) θ is C 1 on ]0, 1[.
ii) θ extends continuously at the origin and let θ(0) := q. (@@)

iii) The majoration θ(x) ≤ q is valid for x ∈]0, ε[ for some ε > 0.

Then define γ := lim sup x→0 + (θ(x) -q).(-Log x).

iv) lim x→0 + θ (x).x.(Log x) = 0.
Note that iii) and iv) are consequences of the condition (@) which can always be assumed when

θ is in C 1 ([0, 1[) or may be written θ(x) = θ 0 (x) -2/(-Log x) with θ 0 ∈ C 1 ([0, 1[).
1.4 The case η(0) = 0 and q = 1

We consider now the standard situation with the conditions q = 1, a = η(0) = 0 and (@). The the computation above in this case gives the inequality:

1/2 xν t -1 .(-Log t) -β-1 .f (t).Z(t).dt ≤ -2 -1 .(Log 2) -β .f (1/2) + β.(Log 2) -β-2 f (1/2) -(Log 2) -β .(θ(1/2) -1).f (1/2) (*bis)
with

Z(t) = β.(β + 1) -Log t -(θ(t) -1).β -θ (t).t.(-Log t) -(-Log t).η(t) (**bis) Now if β.(β + 1) -α 2 > 0 we have Z(t).(-Log t) ≥ β.(β + 1) -α 2 2
when t is small enough, and then we obtain:

1/2 0 t -1 .(-Log t) -β-2 .f (t).dt < +∞ (2.1)
for all β > 0 such that β.(β + 1) > α 2 , so for β > sup{-1/2 + 1/4 + α 2 , 0} when α 2 > -1/4 or for all β > 0 when α 2 ≤ -1/4.

Note that, in this case, it would be enough to assume for θ the conditions i), ii), iii) and iv)-bis lim x→0 + θ (x).x.(Log x) 2 = 0.

in order to obtain (2.1) for β > sup{(γ -1)/2 + (γ -1) 2 /4 + α 2 , 0} when (γ -1) 2 + 4α 2 > 0 or for all β > 0 when (γ -1) 2 + 4α 2 ≤ 0.

1.5 The case η(0) = 0 and q < 1 Assume now that q < 1 and define g(x) := x -s .f (x) where s is a real number that we shall choose later on. We have

f (x) = s.x s-1 .g(x) + x s .g (x) and f (x) = s.(s -1).x s-2 .g(x) + 2s.x s-1 .g (x) + x s .g (x) so x 2 .f (x) + θ(x).x.f (x) -η(x).f (x) =
x s . s(s -1).g(x) + 2s.x.g (x) + x 2 .g (x) + θ(x). s.g(x) + x.g (x) -η(x).g(x)

So the differential inequation (1) for f is equivalent to the differential inequation

x 2 .g (x) + (θ(x) + 2s).x.g (x) -η(x) -s.θ(x) -s.(s -1) .g(x) ≤ 0.

(1 s )

The differential inequality (1 s ) is obtained from (1) by the transformation θ → θ + 2s and η → η -s.θ -s.(s -1).

Now choose s = 1 -q > 0 ; then we obtain for g the differential inequality (1) with θ(x) = θ(x) + 2 -2q, so q := θ(0) = 2 -q > 1 as we assumed q = θ(0) < 1, and with η(x) = η(x) -(1 -q).(θ(x) -q). Note that η is continuous and vanishes at 0; because θ is C 1 we have

α 1 = lim sup x→0 + η(x).(-Log x) = lim sup x→0 + η(x).(-Log x).
Now we apply the case q > 1 to g and we conclude that (2) holds for g and this gives:

1/2 0 t q-2 .(-Log t) -β-1 .t q-1 .f (t).dt < +∞ and then

1/2 0 t -1 .(-Log t) -β-1 .f (t).dt < +∞ (2.2)
for all β > sup{α 1 /(1 -q), 0}.

Note that, in fact, the conditions (@@) on θ are enough to obtain (2.2) for

β > sup{α 1 /(1 -q) + γ, 0}.
1.6 The case η(0) > 0 and q = 1

We consider now the case where a is positive and q := θ(0) = 1. We assume that θ : [0, 1[→ R is still a C 1 and satisfies 2 θ(x) -θ(0) ≤ 0 for x ∈ [0, ε] and we look at positive solutions of the differential inequality

x 2 .f (x) + x.θ(x).f (x) ≤ η(x).f (x) ∀x ∈]0, 1[ (1)
where f is a C 2 function on ]0, 1[. Defining g(x) = x -s .f (x) we have seen that ( 1) is equivalent to

x 2 .g (x) + x. θ(x).g (x) ≤ η(x).g(x) ∀x ∈]0, 1[ with θ(x) = θ(x) + 2s and η(x) = η(x) -s.θ(x) -s.(s -1),
so, choosing for s a solution of the equation s 2 + (q -1).s -a = 0 where q = θ(0), we obtain η(0) = 0, lim sup x→0 + η(x).(-Log x) = α 1 and q = θ(0) = q + 2s.

2 But see the remark at the end of 1.3.

In fact we shall choose the positive root of s 2 + (q -1).s -a = 0: s := s + = -(q -1)/2 + (q -1) 2 /4 + a.

Then q = q + 2s + = 1 + (q -1) 2 + 4a > 1 so we may apply to g the proposition 1.3.1. This gives 1/2 0 t q-2 .(-Log t) -β-1 .g(t).dt < +∞ for all β > 0 such that β > α 1 (q -1), that is to say

1/2 0 t (q-3)/2+ √ (q-1) 2 /4+a .(-Log t) -β-1 .f (t).dt < +∞ for all β > sup{α 1 (q -1) 2 + 4a , 0}.
Note that for a = 0 this computation gives for q > 1 (q -3)/2 + (q -1) 2 /4 + a = q -2 and for q < 1 (q -3)/2 + (q -1) 2 /4 + a = -1.

1.7 The case η(0) > 0 and q = 1

We assume now that the continuous function η satisfies lim sup

x→0 + (η(x) -a).(-Log x) = 0.
where a is a positive real number (so α 1 = 0 and η(0) > 0). Then let the function g defined by the equality f (x) = x s .(-Log x).g(x). Then we have

x 2 .f (x) + θ(x).x.f (x) -η(x).f (x) =
x s .(-Log x). x 2 .g (x) + 2s.x.g (x) -2s.g (x)/(-Log x) + s.(s -1).g(x) + (2s -1).x s .g(x)+ θ(x). x s .(-Log x).x.g (x) + s.x s .(-Log x).g(x) + x s .g(x) -η(x).x s .(-Log x).g(x)

So the inequality x 2 .f (x) + θ(x).x.f (x) ≤ η(x).f (x) is equivalent to the inequality x 2 .g (x)+ θ(x)+2s-2/(-Log x) .x.g (x) ≤ η(x)-s.θ(x)-s.(s-1)- θ(x) + 2s -1 -Log x .g(x)
.

Choosing s := -√ a we obtain as θ(x) := θ(x) + 2s -2/(-Log x) and, if θ satisfies the conditions (@@), the same will be true for θ (but changing γ to γ -2), then

q := θ(0) = 1 -2 √ a and η(x) := η(x) -s.θ(x) -s.(s -1) - θ(x) + 2s -1 -Log x .
Then η satisfies η(0) = 0, as θ(0) = 1 and lim sup

x→0 + η(x).(-Log x) = 2 √ a.
So we are under the hypothesis of the case q < 1 and η(0) = 0 for g, as η has the suitable estimate at 0 + . As α 1 (1 -q) = 1 we obtain that (here γ = -2)

1/2 0 t √ a-1 .(-Log t) -β-1 .f (t).dt < +∞ ∀β > 0 (2.3)

Examples

The following computation shows that the result of proposition 1.3.1 is optimal. Let f (x) := x -λ .(-Log x) β and η(x) := α 1 /(-Log x) + α 2 /(-Log x) 2 for x positive small enough, where α 1 is positive.

Then we have

f (x) = -λ.x -λ-1 .(-Log x) β -β.x -λ-1 .(-Log x) β-1 f (x) = (λ + 1).λ.x -λ-2 .(-Log x) β + β.(2λ + 1).x -λ-2 .(-Log x) β-1 + + (β -1).β.x -λ-2 .(-Log x) β-2
If we define g by the equality

x 2 .f (x) + q.x.f (x) = g(x).f (x) we obtain g(x) := λ.(λ + 1 -q) + β.(2λ + 1 -q).(-Log x) -1 + (β -1).β.(-Log x) -2 .

Assuming that q > 1 and choosing λ + 1 -q = 0 and β.(2λ + 1 -q) = α 1 so λ = q -1 and β.(q -1) = α 1 , we obtain, for q -1 small enough in order that β.(β -1) < α 2 implies g(x) ≤ η(x) for x small enough,

x 2 .f (x) + q.x.f (x) ≤ η(x).f (x)
Remark now that, for the choice β = α 1 /(q -1) the convergence of the integral 1/2 0 x q-2 .(-Log x) -δ-1 .f (x).dx holds if and only if β -δ -1 < -1 that is to say for δ > β = α 1 /(q -1) (compare with proposition 1.3.1).

For q = 1, choosing η(x) := α 2 (Log x) 2 where α 2 is positive, λ = 1 and β := β + , the positive root of (β -1).β = α 2 , the integral 1/2 0

x -1 .(-Log x) -δ-2 .f (x).dx converges if and only if β + -δ -2 > -1 so for δ > β + -1 = -1/2 + 1/4 + α 2 (compare with the conclusion of 1.4).

The following example shows that for q = 1 and η(0) = 0 the differential inequality (1) does not implies the point-wise bound near 0 of a solution by C.(-Log x) for any positive constant C:

For x ∈]0, 1/3[, let f (x) = (-Log x).Log(-Log x). Then f (x) = -1 x .Log(-Log x) - 1 x f (x) = 1 x 2 .Log(-Log x) + 1 x 2 .(-Log x) + 1 x 2
and so, we obtain

x 2 .f (x) + x.f (x) = 1 (-Log x) = η(x).f (x)
where η(x) := 1 (Log x) 2 .Log(-Log x).

.

We have lim x→+∞ (Log x) 2 .η(x) = 0 but f (x)/(-Log x) goes to +∞ when x goes to 0 + .

An application

Our application is devoted to a kind of generalization of Gronwald's lemma for the Laplacian in a punctured ball in R n . We begin by the case n ≥ 3.

Corollary 2.0.1 Let n ≥ 3 be an integer and consider a C 2 function F : B * r → R + on the punctured ball B * r with center 0 and radius r in R n endowed with its natural euclidian norm. Assume that there is a continuous function

η : B r → R with lim sup ||x||→0 η(x).||x|| -2 .(-Log ||x||) = α 1 < +∞ such that F satisfies on B * r ∆F (x) ≤ η(x).F (x) ( 6
)
where ∆ is the standard Laplace operator on R n . Then we have, for any β > sup{α 1 (n -2), 0}

B * r/2
||x|| -2 .(-Log(||x||)) -β-1 .F (x).dx < +∞ (7

)
where dx is the euclidian measure.

Proof. Define the function f :]0, r[→ R + by f (t) := St F (t.σ).dσ where (t, σ) are polar coordinates on R n and where dx = t n-1 .dt ∧ dσ where dσ is the rotation invariant measure on the sphere S 1 with a suitable normalization. Then the function f is C 2 on ]0, r[ and we have

f (t) + n -1 t .f (t) = St ∆(F )(t.σ).dσ ∀t ∈]0, r[ (8)
thanks to the rotation invariance of the Laplacian. Then it is easy to apply the proposition 1.3.1 to the function f with θ ≡ n -1.

The case n = 2 corresponds to the case λ = q = 1. It gives in a similar way the following corollary. where dx is the euclidian measure.

  Let n = 2 and consider a C 2 function F : B * r → R + on the punctured disc D * r with center 0 and radius r in R 2 endowed with its natural euclidian norm. Assume that there is a continuous function η :D r → R with lim sup ||x||→0 η(x).||x|| -2 .(-Log ||x||) 2 = α 2 < +∞ such that F satisfies on D * r ∆F (x) ≤ η(x).F (x) (6)where ∆ is the standard Laplace operator on R n . Then we have, for any β > 0 such that β.(β + 1) > α 2 : * D r/2 ||x|| -2 .(-Log(||x||)) -β-2 .F (x).dx < +∞ (7)

Note that by condition iii) γ ≤ 0.


Florian Schanzenbacher
email: author.florian.schanzenbacher@ratp.fr

Nadir Farhi

Fabien Leurent

Gérard Gabriel

Gérard Gabriel Real-Time

Real-time control of metro train dynamics with minimization of train time-headway variance

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Real-time control of metro train dynamics with minimization of train time-headway variance

Florian Schanzenbächer 1 * , Nadir Farhi 2 , Fabien Leurent 3 , Gérard Gabriel 4

Abstract-We present here a real-time control model for the train dynamics in a linear metro line system. The model describes the train dynamics taking into account average passenger arrival rates on platforms, including control laws for train dwell and run times, based on the feedback of the train dynamics. The model extends a recently developed Max-plus linear traffic model with demand-dependent dwell times and a run time control. The extension permits the elimination of eventual irregularities on the train time-headway. The resulting train dynamics are interpreted as a dynamic programming system of a stochastic optimal control problem of a Markov chain. The train dynamics still admit a stable stationary regime with a unique average growth rate interpreted as the asymptotic average train time-headway. Moreover, beyond the transient regime of the train dynamics, our extension guarantees uniformity in time of the train time-headways at every platform.

I. INTRODUCTION AND LITERATURE REVIEW

This paper is part of a series of works being realized as a research project at RATP, the Paris metro operator. The research project aims to develop a passenger demanddependent traffic management and real-time control algorithm for high density metro lines. The characteristic of the algorithm is its sensibility to passenger demand, i.e. a premier aim is to control the traffic while maintaining a high frequency, satisfying the passenger demand.

The modeling approach is based on former works from [START_REF] Van Breusegem | Traffic Modeling and State Feedback Control for Metro Lines[END_REF], where a traffic control for metro loop lines has been proposed. This approach has served as a basis for the authors of [START_REF] Fernandez | Predicitive traffic regulation for metro loop lines based on quadratic programming[END_REF] which have presented a model predictive control algorithm for metro loop lines aiming for headway regularity one the one side, and schedule adherence, on the other side. Lately, the authors of [START_REF] Schanzenbacher | Train flow optimization: A model predictive and quadratic programming approach (Fluidification du trafic Transilien : approche prédictive et optimisation quadratique)[END_REF] have applied this model predictive control algorithm to a part of a mass transit railway line in the area of Paris. Still on the same basis, the authors of [START_REF] Li | Joint optimal train regulation and passenger flow control strategy for high-frequency metro lines[END_REF], [START_REF] Li | Optimal switched control design for automatic train regulation of metro lines with time-varying passengers arrival flow[END_REF] have proposed a traffic control algorithm for high density metro lines, taking into account the passenger demand.

Starting with a recently developed traffic model for metro lines, with static dwell and run times, in [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF], the authors of [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics in a metro line system with a junction[END_REF] have extended this approach to a metro line with a junction. Applied to a RATP metro line in Paris, this traffic model has precisely reproduced the microscopic characteristics of the line, e.g. the bottleneck segment of the line, the average frequency depending on dwell, run, safe separation times and the number of segments and trains. In [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF], the authors of this paper here, have replaced the static dwell and run times by a dynamic passenger demand-dependent dwell time. A run time control has been introduced, avoiding the amplification of perturbations, i.e. a cascade effect. These traffic models have been shown to be linear in Max-plus algebra, which allows to derive analytic formulas for the average headway and frequency on the line.

In this paper, we refer to a Max-plus Theorem, see [START_REF] Baccelli | Synchronization and linearity: an algebra for discrete event systems[END_REF], [START_REF] Cochet-Terrasson | Numerical computation of spectral elements in maxplus algebra[END_REF], [START_REF] Rob | Railway timetable stability analysis using maxplus system theory[END_REF], having already been used in our preceding works.

In the following, we present an enhanced version of the traffic model which we have proposed in [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF]. It keeps its advantages, i.e. the passenger demand-dependency and the cancellation of a possible cascade effect. As a new element, the here presented model includes a train position control. This means, we modify the model of [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF] such that a minimization of train time-headway variances over the metro line is realized, while maintaining the average frequency on the line, which is a crucial issue for mega cities as Paris, with regard to the very high passenger demand. In the upcoming sections we recall first the modeling approach this paper is based on, before presenting the new traffic model with headway variance minimization. We underpin the practical relevance of this new traffic model with simulation results on a RATP metro line and conclude with an outlook on upcoming work. We consider a linear metro line without junction, as in Fig. 1. The line is discretized in space into a number of segments where the length of every segment is bigger than the train length, as in [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF]. The train dynamics is modeled here by applying run and dwell time control laws under train speed and safe separation constraints. The passenger travel demand is modeled with a static model, considering average arrival rates to every platform. We use the following notations, similar to the ones used in [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF]. For the train dynamics: n the number of segments of the linear metro line. m the number of running trains on the metro line. d k j the k th departure time of trains on segment j. a k j the k th arrival time of trains on segment j. r k j the k th running time on segment j (between nodes j -1 and j). w k j = d k j -a k j the k th dwell time on node j. t k j = r k j + w k j the k th travel time from node j -1 to node j.

II. REVIEW OF THE MODELING APPROACH

g k j = a k j -d k-1 j the k th safe separation time (or close-in time) at node j. h k j = d k j -d k-1 j = g k j + w k j the k th departure time-headway at node j. s k j = g k+bj j -r j . b j ∈ {0, 1}. It is 0 (resp. 1)
if there is no train (resp. one train) on segment j at time zero. bj = 1 -b j .

The average on k and j of the quantities above are denoted r, w, t, g, h and s. We have the following relationships [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF].

g = r + s, (1) t
= r + w, (2)
h = g + w = r + w + s = t + s. (3)
For the passenger demand:

λ ij passenger travel demand from platform i to platform j, when i and j denote platforms, and λ ij = 0 if i or j is not a platform node. λ in i = j λ ij the average passenger arrival rate on origin platform i to any destination platform. λ out j = i λ ij the average passenger arrival rate on any origin platform to destination platform j. α in j average passenger upload rate on platform j. α out j average passenger download rate on platform j.

i λij hi α out j time for passenger download at platform j.

λ in j α in j h j
time for passenger upload at platform j.

We approximate here i λ ij h i /α out j by (λ out j /α out j ) h j . We then define x j as a passenger demand parameter

x j = λ out j α out j + λ in j α in j , (4)
such that x j h j represents the time needed for passenger alighting and boarding at platform j. We notice that capacity limits for trains and platforms are not modeled here. Furthermore, we suppose that the passenger upload and download rate are significantly larger than the travel demand, such that

x j << 1. We also define

X j = x j 1 -x j , (5)
in such a way that x j h j = X j g j . Lower bounds are written r j , w j , t j , g j , s j , h j . The minimum run and dwell time r j , w j are directly given by the infrastructure constraints of the metro system, as well as g j , i.e. the minimum time span between the (k -1) th departure and the k th arrival on segment j. The bounds t j , s j , h j can be calculated based on the lower bounds above, applying relationships [START_REF] Baccelli | Synchronization and linearity: an algebra for discrete event systems[END_REF][START_REF] Van Breusegem | Traffic Modeling and State Feedback Control for Metro Lines[END_REF][START_REF] Cochet-Terrasson | Numerical computation of spectral elements in maxplus algebra[END_REF].

Respectively, upper bounds used here, are denoted wj , ḡj , hj . Let us denote by κ the passenger load capacity per train. The upper bound hj can be fixed as follows.

hj := κ/(

j λ in j ), ∀j, (6)
in such a way that all the passenger arrivals can theoretically be served. Then, based on hj , the upper bounds wj and ḡj are fixed as follows.

wj := x j hj , ∀j, (7) ḡj
:= hj -w j , ∀j. (8)
Let us finally note rj the nominal run time on segment j. rj is given by stability condition [START_REF] Schanzenbacher | Train flow optimization: A model predictive and quadratic programming approach (Fluidification du trafic Transilien : approche prédictive et optimisation quadratique)[END_REF] below and depends on ∆r j . Indeed ∆r j influences the robustness and the frequency of the metro system, i.e. it can either be fixed equal to ∆w j to enhance the average frequency, or it can be chosen greater than ∆w j to reinforce robustness towards perturbations.

B. Review of the Max-plus linear traffic model [11]

We shortly review here the Max-plus linear train dynamics for a metro loop line, with demand-dependent dwell and run time controls [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF]. The dwell and run time control laws are fixed as follows.

w k j = min{x j h k j , wj }, (9)
r k j = max{r j , rj -x j (h k j -h j )}. (10)
The authors of [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF] have shown that under the two following conditions

h 1 j ≤ h1 j = 1/(1 -x j )ḡ j , (11)
∆r j = rj -r j ≥ ∆w j = X j ∆g j , ∀j, (12)
(9) and ( 10) sum to

t k j = r k j + w k j = rj + X j g j . ( 13
)
Then assuming that the k th departure from segment j is realized as soon as the two constraints on the travel and on the safe separation time are satisfied, the train dynamics are written as follows.

d k j = max d k-bj j-1 + rj + X j g j , d k-bj+1 j+1 + s j+1 . (14)
It has been shown in [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF] that the train dynamics ( 14) admit a stationary regime with a unique average growth rate, interpreted as the asymptotic average train time-headway h. Moreover, the latter has been derived analytically by applying basic results of the Max-plus algebra theory [START_REF] Baccelli | Synchronization and linearity: an algebra for discrete event systems[END_REF], [START_REF] Cochet-Terrasson | Numerical computation of spectral elements in maxplus algebra[END_REF], in function of m, the number of running trains on the metro line, X, the passenger demand, and other parameters. We recall this principal result below.

h(m, X) = max                j rj +Xj g j m , max j rj + X j g j + s j , j s j n-m (15)
Under conditions ( 11) and ( 12), the dynamics are Maxplus linear and admits a stationary regime with a unique average train time-headway h(m, X) as shown in (15). These results are very interesting since we guarantee that under the control laws ( 9) and ( 10) on the train dwell and run times respectively, and under the two conditions ( 11) and ( 12), the dwell times respond, i.e. are extended due to the passenger demand, whereas the run times are adapted in order to compensate eventual extensions of the dwell times, to guarantee the overall train dynamics remain stable.

However, we know that the deviation in time of the train time-headway from its average value does not necessarily converge to zero, but it can be periodic around zero. Consequently, irregularities on the headway h k j are not minimized over the k th departures. This is due to the fact that the Maxplus linear system [START_REF] Li | Optimal switched control design for automatic train regulation of metro lines with time-varying passengers arrival flow[END_REF] may admit more than one Max-plus eigenvectors associated to its unique Max-plus eigenvalue.

C. Review of the dynamic programming traffic model [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF], [START_REF] Farhi | Traffic modeling and real-time control for metro lines. Part II -The effect of passengers demand on the traffic phases[END_REF], [START_REF] Farhi | Physical models and control of the train dynamics in a metro line without junction[END_REF] We recall in this section some dynamic models proposed and studied in [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF], [START_REF] Farhi | Traffic modeling and real-time control for metro lines. Part II -The effect of passengers demand on the traffic phases[END_REF], [START_REF] Farhi | Physical models and control of the train dynamics in a metro line without junction[END_REF]. In those models, the train run times are assumed to be constant on all the segments of the metro line. A first studied model fixed the train dwell times as follows.

w k j = max{w j , (λ j /α j )g k j }. (16)
Formula ( 16) models the fact that the train dwell times respect lower bounds w j at every platform, and are eventually extended to embark all the passengers accumulated during the close-in time g k j . It has been recalled in [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF] that the train dynamics with such dwell times is unstable. The model ( 16) justifies the natural instability of the train dynamics, which has already been shown in [START_REF] Van Breusegem | Traffic Modeling and State Feedback Control for Metro Lines[END_REF].

In [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF], [START_REF] Farhi | Traffic modeling and real-time control for metro lines. Part II -The effect of passengers demand on the traffic phases[END_REF], [START_REF] Farhi | Physical models and control of the train dynamics in a metro line without junction[END_REF], a new control law for the train dwell times at platforms has been proposed. The dwell times have the following form.

w k j = max{w j , wj -δ j g k j }, (17)
where δ j are parameters inversely proportional to λ j , in such a way that w k j will be directly proportional to λ j . An important fact here is that w k j is inversely proportional to g k j , i.e. for a long g j , the dwell time at the upcoming platforms is shortened. This is advantageous to minimize headway irregularities without degrading the frequency, as it would be the case for a holding control policy. However, this is contraintuitive from a passenger point of view, since, a retard leads to an accumulation of passengers on the platform, i.e. trains entering the station with a long interval to the preceding train, have to cope with a higher demand.

The authors of [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF] have shown that with dwell time control (17), and under the condition 0 ≤ δ j ≤ 1, ∀j, the train dynamics admits an asymptotic regime with a unique average growth rate interpreted here as the average train time-headway. Therefore the train dynamics are stable. Moreover, they are interpreted as the dynamic programming system of a stochastic optimal control problem of a Markov chain.

III. TRAFFIC MODEL WITH HEADWAY VARIANCE

MINIMIZATION

On the one side the Max-plus linear traffic model of section II-B from [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF] guarantees stability and controls both dwell and run times to take into account the passenger travel demand. However, it does not harmonize train departure time intervals on the line.

On the other side, the dynamic programming model of section II-C from [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF], [START_REF] Farhi | Traffic modeling and real-time control for metro lines. Part II -The effect of passengers demand on the traffic phases[END_REF], [START_REF] Farhi | Physical models and control of the train dynamics in a metro line without junction[END_REF] guarantees stability but controls only the dwell times at platforms without controlling the run times. However, this model is better in term of harmonization of train departure time intervals, comparing to the Max-plus one.

The new model we propose here extends the Max-plus model ( 14) of section II-B to a dynamic programming model as the model of section II-C. Consequently, our new model benefits from the advantages of both models.

• Accounting for passenger demand with (partially) demand-dependent dwell times and run time control, • Minimizing train time-headway variances on the platforms via a dwell time control without degrading the average frequency. We consider here the same train run time control law [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics in a metro line system with a junction[END_REF], and propose a new train dwell time control law which modifies (9) as follows.

w k j = min (1 -γ j )x j h k j , wj , (18)
with 0 ≤ γ j ≤ 1.

We notice that by activating the control (0 < γ j ≤ 1) in case of a perturbation on the train headway, we limit excessively long dwell times, which would have been a direct consequence of a long headway.

Let us use the following additional notations.

• ∆w j := wj -w j .

• ∆g j := ḡj -g j .

• ∆r j := rj -r j . It is then easy to check that ∆w j = 1/(1 -x j )∆g j .

Proposition 1: If for all j, h 1 j ≤ hj and ∆r j ≥ ∆w j , then t k j = rj + X j g j -γ j X j g k j , ∀j. Proof: The proof is by induction. It is similar to the one of Theorem 1 in [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF].

The train dynamics is then written as follows.

d k j = max      (1 -δ j )d k-bj j-1 + δ j d k-1 j + (1 -δ j )(r + X j g j ), d k-bj+1 j+1 + s j . (19
) with δ j := (γ j x j )/(1 + γ j x j ).

Let us notice that the dynamics (19) extends [START_REF] Li | Optimal switched control design for automatic train regulation of metro lines with time-varying passengers arrival flow[END_REF], since in the case where δ j = 0 (i.e. γ j = 0), (19) coincides with [START_REF] Li | Optimal switched control design for automatic train regulation of metro lines with time-varying passengers arrival flow[END_REF]. This is a direct consequence of the fact that the train dwell time control law (18) extends [START_REF] Rob | Railway timetable stability analysis using maxplus system theory[END_REF].

It is easy to see that if m = 0 (zero trains) or if m = n (the metro line is full of trains), then the dynamic system ( 19) is fully implicit. Indeed, if m = 0, then b j = 0, ∀j, i.e. the first term of the maximum operator in ( 19) is implicit for every j. Similarly, if m = n, then bj = 0, ∀j, i.e. the second term of the maximum operator in ( 19) is implicit for every j. We notice that in both cases, m = 0 and m = n, no train movement is possible. On the other side it is not difficult to check that if 0 < m < n, then the dynamic system is implicit but triangular, i.e. there exists an order on j of updating the variables d k j in such a way that the system will be explicit. In fact, this order corresponds to the one of the train movements on the metro line. In the following, we consider only the case 0 < m < n. Therefore, the dynamic system (19) admits an equivalent triangular system.

The train dynamics (19) can be written as follows.

d k j = max u∈U M u d k-1 j + N u d k j + c u j , ∀j, k, (20)
where M u , u ∈ U and N u , u ∈ U are two families of square matrices, and c u , u ∈ U is a family of column vectors. Moreover, since 0 ≤ δ j < 1, ∀j by definition, then M u ij ≥ 0, ∀u, i, j and N u ij ≥ 0, ∀u, i, j. Furthermore, we have j M u ij + N u ij = 1, ∀u, i. The equivalent triangular system of system (19) can be written as follows.

d k j = max u∈U M u d k-1 j + cu j , ∀j, k, (21)
where M u , u ∈ U is a family of square matrices, and cu , u ∈ U is a family of column vectors, which can be derived from M u , N u and c u , u ∈ U. Moreover, we still guarantee M u ij ≥ 0, ∀u, i, j and j M u ij = 1, ∀u, i. By consequent, the system (21) can be seen as a dynamic programming system of an optimal control problem of a Markov chain, whose transition matrices are M u , u ∈ U associated to every control action u ∈ U, and whose associated rewards are cu , u ∈ U.

Theorem 1: If 0 < m < n, then the dynamics (19) admits a stationary regime, with a unique asymptotic average growth vector (independent of the initial state vector d 0 ) whose components are all equal to h, which is interpreted here as the asymptotic average train time-headway.

Proof: We give only a sketch of the proof here. Since the dynamics (19) are interpreted as the dynamic programming system of a stochastic optimal control problem of a Markov chain, it is sufficient to prove that the Markov chain in question, which is acyclic since 0 < m < n, is irreducible for every control strategy. This is not obvious from the dynamics (19). We need to show it on the equivalent triangular system (21). An alternative proof is the one of Theorem 5.1 in [START_REF] Farhi | Traffic modeling and real-time control for metro lines[END_REF].

Theorem 1 does not give an analytic formula for the asymptotic average train time-headway h. However, it guarantees its existence and its uniqueness. Therefore, h can be approximated by numerical simulation based on the value iteration, as follows.

h ≈ d K j /K, ∀j, for a large K. (22)
Let us now go back to the objective of the extended model we propose here, i.e. our model (19) extending the model ( 14) of [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF]. We give here the main ideas explaining why the control we proposed permits to minimize the asymptotic average train time-headway variance. Rigorous arguments with proved theorems will be given in our forthcoming articles.

As mentioned above, if γ j = 0, ∀j, the dynamics (19) are Max-plus linear, and are equivalent to [START_REF] Li | Optimal switched control design for automatic train regulation of metro lines with time-varying passengers arrival flow[END_REF]. In this case, if we write the dynamics (19) or equivalently [START_REF] Li | Optimal switched control design for automatic train regulation of metro lines with time-varying passengers arrival flow[END_REF] under the triangular form (21), then the associated matrices M u are boolean circulant matrices. The latter may then have eigenvalue 1 with multiplicity bigger than 1. Therefore, although the average growth rate of the dynamics is unique and independent of the initial state d 0 , the asymptotic state d k (up to an additive constant) may depend on the initial state d 0 .

However, in case of the dynamics (19) with 0 < γ j ≤ 1 for some j, the application of γ j will force the dynamics to converge a stationary regime where the activated matrices M u are the ones having eigenvalue 1 as a simple one, i.e. with multiplicity 1. In this case, d k will converge, up to an additive constant, independent of the initial state d 0 . Moreover, the asymptotic state d will correspond to the case where the train time-headway is the same at all the platforms.

We will show in the next section by means of numerical simulation, that in this case, the train dynamics converge to a state of traffic where the train departures are harmonized. In other words, while the train time-headway converges to its asymptotic value h, the variance in time of the train time headways, i.e. its deviation with respect to its average value h, converges to zero.

With regard to dwell time equation ( 18), we notice that in case γ j > 0 where the control is applied, we do not only ensure train headway convergence at platforms, but are also sure not to degrade the train frequency. Indeed, with this, dwell times are slightly shortened to allow minimization of variances on the train time-headway. This means, in case of a big interval, train dwell times are no longer extended excessively, but are expanded just a little bit. On networks with a very high demand, as for the RATP network in Paris, this is highly advantageous to a holding control, because it allows a passenger load harmonization while maintaining a high train frequency.

IV. SIMULATION RESULTS

In Tab. I we present some simulation results for different values of parameter γ j . First, we depict the results for the Max-plus dynamics, i.e. γ j = 0. Second, we fix γ j = 0.1, which is close to the Max-plus dynamics. Third, γ j degrades linearly from 0.5 to 0, over the simulation horizon K. The simulation horizon K is fixed accordingly to the time span, after which the train positions should be harmonized. It has been chosen to K = 80 departures here. In the left column, we show the train trajectories, i.e. the departure times over all the segments on the metro loop line. In the right column, we give the last observed headways h(γ j ) = d(γ j ) K j -d(γ j ) K-1 j over all the segments j on the line, indeed a good way to analyze the efficiency of the proposed control strategy, aiming to minimize headway variance. The first row presents the results for γ j = 0, i.e. without control where the dynamics (19) are a Max-plus linear system, as in [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF]. It can easily be seen, that, depending on the initial train positions, we observe some areas with more trains on the line, and others with less, i.e. the headways are irregular. Moreover, this irregularity persists over the simulation. In fact, it has been shown in [START_REF] Schanzenbacher | A discrete event traffic model explaining the traffic phases of the train dynamics on a linear metro line with demand-dependent control[END_REF] that the train time-headways along the line do not converge. They rather reach a periodic regime which makes their asymptotic average convergent. We notice that the train time-headways do not diverge although dwell times are modeled depending on the headway, which remain constant with the applied run time control. The graph giving the last observed headway shows clearly that important irregularities on the headway persist without control, depending on the number of trains, ranging from less than 5 minutes up to over 15 minutes.

The second row depicts the results for a static control where γ j = 0.1 for all j which include a platform. In the left graph we see some irregularities at the beginning of the simulation, i.e. at the bottom, which tend to be minimized over time. In fact, with regard to the final headways in the right graph, it can be seen that headways clearly converge along the line. However, due to the low γ j = 0.1, larger initial variances cannot be minimized fast enough, i.e. they still can be observed at the end of the simulation. Variances finally range between less than 5 minutes and 10 minutes, i.e. a maximum variance of around 5 minutes.

In the last row, we give the results for a dynamic linear control where γ j = 0.5 -(0.5/K)k for all j with a platform. The graph in the left column shows that some headway irregularities at the beginning are quickly being minimized and headways converge. In the right column, the graph giving the final headways observed along the line, shows that headways converge more quickly than for the static control. For the same simulation horizon as in the static case, the final headway variance does not exceed 2.5 minutes, i.e. headways range between slightly less than 5 minutes and around 7.5 minutes. Finally, the dynamic control is advantageous because it applies stronger controls right after a perturbation in the transient regime, but returns to a γ j = 0 once the perturbation has been absorbed, i.e. once the stationary regime has been reached. From a passenger point of view, this is preferable since for a γ j = 0, dwell times respond once more fully to the time needed for passenger de-and embarking.

V. CONCLUSION AND FUTURE WORKS

Our series of traffic models for metro lines, ranging from a standard model with static run and dwell times, over an enhanced version with demand-dependent dwell times and a run time control, has been completed here by a realtime control with minimization of the train time-headway variance. The simulation results have underpinned the effectiveness of the here proposed traffic control. The figures clearly show headway convergence for the here proposed control. Moreover, they strongly suggest that a dynamic control, where a stronger control is applied right after a perturbation, followed by a fade-out time, leads to better results, i.e. faster headway variance minimization than a light but permanent control.

To complete this control model series, this paper will be followed by an application of this control policy to the traffic model for a metro line with a junction. At this point, we assume to have a complete toolbox of traffic models which will allow an on-site implementation on the RATP network in Paris. Possible applications are metro loop lines as well as metro lines with junctions, where the here presented model ensures headway regularity while maintaining a desired train frequency, with passenger demand-dependent dwell times and run time control. There is always room for further improvement, e.g. for a passenger stock model in the trains and on the platforms, as well as an extension towards more complicated operational concepts, i.e. including intermediate terminus stations and non-stopping policy at selected stations.

Fig. 1 .

 1 Fig. 1. Schema of a metro loop line and the corresponding notation.

  arXiv:1809.02876v1 [math.OC] 8 Sep 2018 A. Notations