Serge Mora
email: smora@univ-montp2.fr

Ty Phou

Jean-Marc Fromental

Yves Pomeau

Gravity driven instability in elastic solids

Keywords: numbers: 46.32.+x, 46.25.-y, 47.20.Ma, 83.80.Kn

We demonstrate the instability of the free surface of a soft elastic solid facing downwards. Experiments are carried out using a gel of constant density ρ, shear modulus µ, put in a rigid cylindrical dish of depth h. When turned upside down, the free surface of the gel undergoes a normal outgoing acceleration g. It remains perfectly flat for ρgh/µ < α * with α * 6, whereas a steady pattern spontaneously appears in the opposite case. This phenomenon results from the interplay between the gravitational energy and the elastic energy of deformation, which reduces the Rayleigh waves celerity and vanishes it at the threshold.

Many materials such as biological tissues can withstand huge elastic deformations of more than several hundred percent. The amplitude of the stress is then of the order of the elastic modulus, a situation commonly encountered with soft materials. Specific and fascinating patterns, reminiscent of those that can be seen in hydrodynamics, can then occur spontaneously [1][2][3][4][5][6][7][8][9][10].

Since both soft elastic solids and liquids are capable of undergoing large deformations, and are often subjected to forces with a common origin, eg capillary forces [6,11], it is likely that some mechanical instabilities can be shared, to a certain extent, by these two kinds of continuous media [3]. Out of the many instabilities experienced by liquids, the Rayleigh-Taylor instability (RTI) [12][13][14] is outstanding because it is easy to understand, not too difficult to rationalize and also important in many technological and physical situations. The dispersion relation for regular gravity waves on a deep ocean reads ω = √ gk, where g is the downward gravity acceleration, ω/(2π) is the wave frequency and k its horizontal wave number. As often noticed, if one turns the gravity upward, that is if one changes the sign of g, ω becomes purely imaginary ±i √ -gk,

showing the existence of fluctuations growing exponentially with time. These fluctuations do not saturate and yield ultimately fingers of liquids in free fall. If one considers, as we do below, a soft solid in air with its surface turned downward, there are a priori good reasons to believe that some sort of RTI will set in. To figure it, consider a horizontal elastic slab of thickness h subjected to the gravity of Earth g, the upper surface being fixed to a rigid body, the lower one being free (Fig. 1).

A sinusoidal perturbation ζ = ε sin(kx) of the surface height (with x an in-plane coordinate) causes a reduction in the gravitational energy per unit area equal to 1 2λ λ 0 ρgζ 2 dx = 1 4 ρgε 2 , where ρ is the mass density of the elastic medium and λ = 2π/k is the wave-length. The corresponding elastic energy cost per unit volume scales as the shear modulus µ times the mean squared strain.

In the long wave limit (kh 1) the strain scales as /h: the sample is vertically squeezed from length h to hε above a trough of the wave (region (a) of Fig. 1), it is vertically stretched above a peak from h to h + ε (region (b)), and the deformation varies progressively in between (region (c)). Finally, the mean elastic energy per unit volume scales as h • µ(ε/h) 2 . Comparing the two contributions of the total energy, it appears that the Rayleigh-Taylor buoyancy overcomes elasticity beyond an instability threshold ρgh/µ = α * where α * , the dimensionless proportionality constant for the elastic energy (with a factor of four), is to be found (See Supplement-B for the complete calculation).

In the same situation a thin layer of liquid is always unstable which is equivalent to set to zero the shear modulus in the previous estimate. Therefore, contrary to liquids RTI in a solid has a well defined threshold for layers of finite thickness. Beyond it, the deformation increases up to a finite value for which the elastic cost balances the buoyancy gain: a steady state of equilibrium is then reached.

Although RTI in solids is expected to play a role in many fields such as biology, geology [15,16] or astrophysics [17], both a direct observation and a clear characterization are missing. In [18,19] and [20], a flat metal plate whose thickness is initially periodically modulated with a low amplitude is accelerated by expanding detonation products. The growth of the initial perturbation is observed through the use of x-ray shadow-graphs. It was found to be governed by the yield strength of the elasto-plastic material, the initial amplitude and the plate thickness. More recently, yogurts with a sinusoidal perturbation at the surface were put in a mold and accelerated using a linear electric motor. The stability regions of this elastic-plastic material have been investigated in terms of acceleration, amplitude and wavelength of the initial perturbation [21]. In both cases (flat metal plate and yogurt), the observations consist in evolving states which are clearly associated to plastic deformations of pre-existing periodic ripples at the free surface. Schematically, when the acceleration exerts a strong enough stress on the ripples, the yield stress of the material is overcome at the ripples extremities which begin to flow. In the experiments of [21], the case of an initially flat surface has been briefly investigated. The authors reported the existence of a non stationary surface instability and related its nucleation to the elastic (reversible) deformations of the material.

Their conclusion seems erroneous insofar it is based on a comparison with a theoretical expression valid for samples whose height is much higher than the wavelength, which is not the case in their experiments. The results obtained in the present paper demonstrates unambiguously that if the phenomenon observed by these authors were a consequence of RTI for an elastic solid, their observations would have been different. It is therefore likely that the observed phenomenon is a consequence of the plastic properties of the investigated material.

Following the pioneering analytic work of Drucker [22], RTI in plastic solids has been modeled in the visco-elasto-plastic approximation in order to simulate the growth in amplitude of initial sinusoidal perturbations [23][24][25]. RTI for purely elastic plates with an initially flat surface has been analytically studied by Plohr and Sharp [26], whose results have been generalized few years after [27,28]. They predict for each value of the acceleration the existence of a critical perturbation wave-length beyond which the flat surface is unstable. As a consequence, an elastic plate is always unstable, provided its dimension are large enough compared to the unstable wavelengths. This is in contrast with the instability studied in this letter. On the other hand, Bakhrakh and Kovalev [START_REF] Bakhrakh | Material 5 All-Union[END_REF] have analytically studied the case of an accelerated elastic half space and found that it is unstable with respect to any perturbations with a wavelength larger than 4πµ/(ρg).

We report below the experimental observation of an instability occurring on the surface of a heavy ideal elastic solid pointing downwards. This instability occurs above a threshold and results in steady patterns. At threshold, elasticity exactly counterbalance buoyancy for an infinitesimal perturbation of the free surface. This phenomenon is closely related to Rayleigh waves [START_REF] Rayleigh | [END_REF] since the phase velocity of elastic surface waves decreases as the inwardly gravity increases. This provides a physical interpretation of the instability we have demonstrated, insofar it occurs when the gravity is strong enough to make fully vanish the phase velocity. This viewpoint leads us straightforwardly to calculate the growth rate of the instability.

In our experiments, we use aqueous polyacrylamide gels consisting in a loose permanent polymer network immersed in water. The density of this incompressible elastic material is almost equal to that of water. It behaves as an elastic solid for strains up to several hundreds of percent (Supplement-A). The shear modulus can be tuned over a wide range by varying the concentrations in monomers and cross-linkers, or the temperature. In our experiments, it lies between 30 and 150 Pa. It is measured through indentation tests (Supplement-A). The reagents generating the gel are dissolved in ultra-pure water and poured into the brim in a cylindrical dish whose walls are covered with a thin layer of Velcro loops to prevent any further detachment. After the gel is made and its shear modulus measured, the dish is flipped upside down. Various methods have been tested: (i) reversal when the system is immersed in water (density close to that of the gel), the container is then gently removed out of the water keeping horizontal the free surface; (ii) reversal carried out in air but with a rigid plate keeping flat the surface during inversion. The plate is then gently removed; (iii) direct and fast flipping of the system in air without any special care. The three methods lead to identical results. The surface of the thinnest or the hardest samples remains perfectly flat (Fig. 2-a). In a narrow range of shear moduli and heights non propagating undulations grow spontaneously at the free surface of the gel and remain permanently (Fig. c,d). For lower shear moduli or for greater thicknesses, several cuvettes appear next to each other at the surface, and remain permanently. For a constant thickness, their number (from one to seven in our experiments) and their size depend on the shear modulus (Fig. f,g). In any cases, flipping again the container (so that the free surface is horizontal and upward) leads to the perfectly flat surface we started from. In addition, successive reversals lead to the same observations, except for a particular sample for which the number of cuvettes is either four, either seven (Fig. 2-f,g). We infer that the shear modulus of this sample corresponds to a threshold for the number of cuvettes so that the final configuration of the system is driven by uncontrolled external disturbances. This point is discussed at the end.

To obtain quantitative information about the surface deformation, a regular light grid is projected about the free surface (Fig. 3-a). The gel being transparent and the bottom of the container being white, the observed image of the grid results from one refraction followed by one reflection and another refraction. If the free surface is flat, this image corresponds to the grid without geometric distortion (Fig. 2-a). It is warped if the free surface is deformed, this distortion is the bigger the surface is more deformed (Fig. 2-b-g). To measure the distortion, a rectangular lattice is fitted with the recorded images using the least squares method. The fitting parameters are a translation, the parameters and the orientation of the lattice, and a possible quadratic distortion taking into account the (small) radial decentering optical distortion (Fig. 3-b). Fig. 3-c shows this deviation plotted as a function of µ for samples having the same thickness h = 2.75 cm. Fitting functions a + b(µµ * ) c with these data gives µ * = 44.6 ± 1.8 Pa, demonstrating the existence of an instability threshold at µ * (Fig. 3-c). This threshold corresponds to a critical dimensionless acceleration α * = ρgh µ * = 6.05 ± 0.25. We expect that the onset of instability will show up when the frequency of a mode of propaga-tion of elastic waves [START_REF] Rayleigh | [END_REF][31][32] at finite wavelength becomes zero. We consider an infinite layer of a heavy and incompressible elastic medium of thickness h with a free downwardly facing surface with air, the other surface being fixed on a rigid substrate. We also consider a plane wave propagating in the in-plane x-direction (see Fig. 3) with the (small) displacement u = u(y)e iωt-kx (boldface being for vectors). Putting this displacement field in the equations of motion for an isotropic and incompressible heavy elastic medium [START_REF]Theory of elasticity[END_REF][START_REF] Kuipers | [END_REF] with the boundary conditions described just above, we obtain after linearization a condition for the dimensionless frequency ω = ωh µ/ρ and the dimensionless wave number k = kh (see Supplement-C):

det        0 2 k2 0 s2 + k2 -ω 2 + 2 k2 -α k 2 ks -α k k cosh k -k sinh k s cosh s -s sinh s -k sinh k k cosh k -k sinh s k cosh s        = 0, (1)
for s2 = k2 -ω2 > 0. In Fig. 4-left, ω is plotted from Eq.1 as a function of k for various values of α = ρgµ h . The curves exhibit a local minimum for α > 4.5, resulting in two possible wavelengths for one frequency. The frequency at the minimum becomes zero for α = α * = 6.223 • • • : the propagation speed of the waves is then zero and a sinusoidal perturbation of the surface with the corresponding wave number is stationary. The surface is then linearly unstable. The theoretical value of α * is in good agreement with experimental observations (Fig. 3). The finite size of our samples therefore has no significant effect on the threshold value, nor the gel-air interfacial tension, which is consistent with the calculation of Supplement-B. Furthermore k = 2.12 with h = 2.75 cm corresponds to 8 cm for the wavelength, consistent with snapshots (b) and (c) of Fig. 2. For α > α * we find ω 2 < 0. Writing ω = iΩ, we obtain the growth rate Ω of the instability (Fig. 4right). We find Ω The preceding theoretical study applies for infinitesimally small strains, ie near the threshold.

Well beyond it, the final patterns are obtained after a substantially longer duration (a few seconds for patterns of Fig. 2-e,f,g). A surprising and striking non linear feature of the instability is the difference in the observed patterns between snapshots (f), (g) and (h) of Fig. 2, all the three corresponding to the same shear modulus with the same container size. The first two are directly as introducing a kind of free choice in the evolution of a system shows up here. It implies that the ultimate state reached after such an instability depends not only on the growth of the unstable structure itself but also on uncontrollable or at least hard to control small effects, like various inhomogeneities in space and time. This is clearly evidenced in our experiment.

We have shown that RTI exists in real elastic solids and that it can be observed in everyday's gravity field in soft hydrogels. The instability threshold depends on the shear modulus, the thickness and the density of the sample. Our experimental set-up with soft elastic gels has enabled a quantitative comparison with a linear theory. This has allowed us to identify the basic ingredients of this instability and the way it appears. Measuring the dispersion relation of surface waves can be a mean to detect the proximity of the threshold, and therefore to predict an impending change.

These results open the way for further fundamental studies, for instance concerning the dynamic formation and the large-scale organization of the patterns, which are both of great importance for non linear physics and morphogenesis. RTI in solids should also be found in more complex situations, such as biology, geology and industrial processing, with visco-plastic, visco-elastic or non-isotropic materials. Moreover, the instability is expected to occur in more extreme conditions (high accelerations, strong and non-uniform gravitational fields) where the direct observation is hardly possible. We believe that our work lays foundation to address such more complex cases.

The authors are indebted to E. Bouchaud and M. Destrade for interesting discussions. This work has been supported by ANR under Contract no. ANR-2010-BLAN-0402-1 (F2F). face being here and thereafter for vectors). In this Lagrangian framework, the actual location in space is parametrized by the coordinates of the preimage in the undisturbed (rest) state.

Incompressibility is imposed by writing that the determinant of the first derivatives of R(x, y) is equal to one, namely that

X ,x Y ,y -X ,y Y ,x = 1, (2)
where X ,x is for ∂X ∂x , etc. Let Y = ζ(X) be the equation of the surface in the disturbed state. The boundary condition (b.c.) on the surface is the condition that the pair that a quantity is evaluated on the free surface (y = 0, Fig. 6).

(X ∂ (x) = X(x, 0), Y ∂ (x) = Y (x, 0)) is a parametric representation of the curve of equation Y ∂ = ζ(X ∂ ),
Assuming a capillary energy proportional to γ, the energy of the system reads:

E = dxdy (W + ρgy) + γ dx X 2 ∂,x + Y 2 ∂,x , (3)
with W the strain energy density function of the elastic material [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF]. Imposing zero variation of E with respect to X(.) and Y (.) yields the equilibrium equations, both in the bulk (Cauchy-Poisson equations) and at the surface (boundary conditions). The general equations have been previously [39] established for the neo-Hookean elasticity, known to describe fairly well the soft solids, or hydrogels, used for experiments. It is the natural extension of Hooke's law and it does not introduce any new physical parameter. From [39], the Cauchy-Poisson equations reads

∂ ∂x (µX ,x -qY ,y ) + ∂ ∂y (µX ,y + qY ,x ) = 0, (4)
∂ ∂x (µY ,x + qX ,y ) + ∂ ∂y (µY ,y -qX ,x ) = ρg, (5)
and the boundary conditions at y = 0 are:

µX ,y + qY ,x -γ ∂ ∂x   X ∂,x X 2 ∂,x + Y 2 ∂,x   = 0, (6)
µY ,y -qX ,x + γ ∂ ∂x   Y ∂,x X 2 ∂,x + Y 2 ∂,x   = 0, (7)
where the Lagrange multiplier q(x, y) allows to impose the condition X ,x Y ,y -X ,y Y ,x = 1 everywhere.

We consider the infinitesimal perturbation of the base state

         X(x, y) = x + εX (y)e ikx
Y (x, y) = y + εY(y)e ikx q(x, y) = q (0) + εq(y)e ikx with ε 1. From Eq.4-5, q (0) = Cρgy (where C is a constant of integration) and Eq.7 imposes

C = µ.
From incompressibility (Eq.2):

ikX + Y ,y = 0. (8)
The Cauchy-Poisson equations (Eq.4-5) read:

-

k 2 µX -ikq + µX ,yy -ikρgY = 0 ⇒ Y ,yyy k 2 -Y ,y - ρg µ Y = q µ , (9)
-k 2 µY + µY ,yy -ikρgXq ,y = 0 ⇒ k 2 Y -Y ,yy + ρg µ Y ,y + q ,y µ = 0.

Substituting the Lagrange multiplier, we obtain:

Y ,yyyy -2k 2 Y ,yy + k 4 Y = 0. (10)
The boundary conditions y = 0 (Eqs.6-7) simplify:

X ,y + Y ,x = 0, Y ,y -X ,x - q µ = - γ µ Y ,xx ,
and with Eq.8 one gets for y = 0:

Y ,yy + k 2 Y = 0, (11)
2Y ,y - q µ = γ µ k 2 Y. (12)
From Eqs.9 and 12 one obtains

3Y ,y - 1 k 2 Y ,yyy = - ρg µ + γ µ k 2 Y. (13)
A zero horizontal displacement at y = h is imposed because the elastic medium is assumed to be bonded to the substrate. Hence the boundary conditions at y = h are Y(h) = 0 and X (h) = 0, or equivalently (from Eq.8)

Y(h) = 0 and Y ,y (h) = 0. ( 14
)
The solutions of Eq.10 are of the form Y = (a+by)e -ky +(c+dy)e ky . Inserting this expression in Eqs.11-13 (for y = 0) and Eqs.14 (for y = h) one obtains: This system has a nonzero solution (a, b, c, d) for:

                 ka -b + kc + d = 0 ρg µ - γ µ 2 -2k a + ρg µ - γ µ k 2 + 2k c =
1 -1 1
ρgh µ + 1 + 2 k2 -2 + ρgh µ 1 k e -2 k = 0, (15)
or equivalently, in

α = 2 cosh(2 k) + 2(1 + 2 k2 ) + γ µh k(sinh(2 k) -2 k) sinh(2 k) k -2 (16)
with k = kh and α = ρgh µ . The flat interface is neutrally stable against an infinitesimal perturbation of dimensionless wave number k for a dimensionless acceleration α = ρgh µ given by the condition of Eq.16. For α larger than the value calculated from Eq.16 the flat interface is unstable against an infinitesimal perturbation of dimensionless wave number k.

Let us first suppose that γ/µ h. In this case The minimum of α depends on the surface tension: its value increases as γ µh is larger, and the corresponding value of k is smaller (Fig. 7). As expected, the surface tension stabilizes the interface, an effect which is more pronounced as the wavelengths are shorter since the most unstable wave number is shifted to the smaller values when the surface tension is increased.

α = 2 k 2 k2 + cosh(2 k) + 1 sinh(2 k) -2 k . ( 17
The polyacrylamide gel-air interfacial tension is of the order of magnitude of the water-air surface tension as the gel mainly consists of water. Therefore, γ ∼ 70 mN/m in our experiments.

The sample thickness being h ∼ 3 cm and the measured modulus at the threshold µ ∼ 45 P a, one obtains γ µh ∼ 0.05. From Eq.16, the shift of α at the threshold is not expected to be significant (see Fig. 7). Therefore surface tension effects are not expected to play a role in our experiments at theorem of Helmholtz we write :

u x = ∂φ ∂x + ∂ψ ∂y (18)
u y = ∂φ ∂y - ∂ψ ∂x (19)
where φ and ψ are functions of x, y and t. φ and are solutions of equations (see Eq.3.3 and Eq.3.5.1 of [START_REF] Kuipers | [END_REF]):

∆φ = 0 and ∆ψ = 1 c 2 2 ∂ 2 ψ ∂t 2 (20)
with c 2 = µ/ρ. Considering a sinusoidal wave of frequency ω/2π and wavelength k/2π propagating in the x-direction, φ and ψ can be chosen as (see Eq.3.9 in [START_REF] Kuipers | [END_REF]):

φ = (A 1 cosh ky + A 2 sinh ky) e i(ωt-kx) (21)
= (B 1 sinh sy + B 2 cosh sy) e i(ωt-kx)

The boundary conditions at the free surface y = 0 are (Eqs.3.7 in [START_REF] Kuipers | [END_REF]):

   2 ∂ 2 φ ∂x∂y + ∂ 2 ψ ∂y 2 -∂ 2 ψ ∂x 2 = 0 ρg ∂φ ∂y -∂ψ ∂x + ρ ∂ 2 φ ∂t 2 + 2µ ∂ 2 φ ∂y 2 -∂ 2 ψ ∂x∂y = 0 (23)
A zero vertical displacement at y = h: u y (h) = 0 is imposed. But contrary to [START_REF] Kuipers | [END_REF] where the elastic material could freely slide parallel to the substrate, here a zero horizontal displacement at y = h is also imposed because the elastic medium is assumed to be bonded to the substrate. Hence Following [START_REF] Kuipers | [END_REF] we introduce the dimensionless parameters:

k = kh, ω = h c 2 ω, α = ρgh µ , s = hs
The condition for the linear system (Eq.25) has nonzero solutions, ie for a wave with a wave number k and a frequency ω/2π can propagate, is: From Eq.20, s 2 = k 2 -(ω/c 2 ) 2 . Hence:

0 -2i k 0 s2 + k2 2
s = k2 -ω2 if k > ω s = i ω2 -k2 = is if k < ω
If k > ω the propagation condition writes: Note that taking ω = 0 in Eq.26 (and hence s = k) yields Eq.17.

FIG. 1 :

 1 FIG.1: Scheme of a sinusoidal disturbance of a downwardly facing and initially flat surface of a heavy elastic slab. The other surface is fixed on the rigid substrate. The energy change due to this disturbance can be positive or negative depending on the value of the dimensionless ratio ρgh/µ with ρ the mass density, g the (gravitational) acceleration, h the thickness and µ the elastic modulus.

FIG. 2 :

 2 FIG. 2: Views of the downwardly facing free surface of gels with different shear moduli. (a) µ = 78± 0.5 Pa; (b) 44 ± 0.5 Pa; (c) 43.3 ± 0.5 Pa; (d) 42.8 ± 0.5 Pa; (e) 41.0 ± 0.5 Pa; (f,g,h) 40.0 ± 0.5 Pa. The cylindrical dish is 18 cm diameter and 2.75 cm deep. (a-g) Steady patterns obtained just after reversal. (h)The sample temperature was 50 degrees when reversed. The snapshot is taken after the sample has cooled to room temperature.

FIG. 3 :

 3 FIG.3: Quantitative analysis of the surface distortion. (a) Experimental setup. The distance mirrorprojector-sample is 1.5 m while the dish diameter is 18 cm. (b) Full circles: intersections of the distorted lines of the grid observed on a sample (µ = 41.5 ± 0.5 Pa). Solid lines defined the grid that best fits the observed one. (c) Square root of the mean squared error with the grid that best fits the observed one as a function of µ with h = 2.75 cm (full circles). The solid line is the best fit with the power law a+b(µ * -µ) c . We find µ * = 44.6 ± 1.8 Pa.

1 .8µ ρh 2 √

 12 αα * for α ≤ 1.3α * , corresponding to a characteristic time (1/Ω) ranging from 0.2 (sample (b) of Fig.2) to 0.1 s (sample (g)). Unfortunately, such a characteristic time cannot be experimentally measured since it is shorter than the duration required to place the sample.

1 FIG. 4 :

 14 FIG.4: Predictions for waves propagating at the surface of a heavy elastic material. Left: dimensionless frequency squared ω 2 h 2 µ/ρ as a function of the dimensionless wave number kh for different values of α.For α < α * , the flat surface is stable and ω is the Rayleigh frequency. For α > α * the flat surface is unstable and the growth rate of the most unstable mode is Ω = -ω 2 min , ω min being the minimum of ω. Right:Dimensionless growth rate of the mode k with the maximum growth rate obtained from Eq.1 as a function of the gap with the instability threshold (solid line).

FIG. 5 :

 5 FIG. 5: (a) Loss and storage moduli. The amplitude of the sinusoidal strain (frequency 2π/ω) is 1%. (b) Stress versus shear strain. (c) Indentation tests using a flat-ended cylindrical punch, 3 mm diameter. Empty circles : Total normal force applied to a cylindrical indenter as a function of the indentation depth at the surface of a gel placed in a cylindrical dish facing upwardly and immersed in water. Dashed line : Buoyancy acting on the indenter calculated from the immersed volume. Full circles : Elastic contribution obtained by subtracting buoyancy to the total force. Solid line: linear fit for the elastic contribution leading to a shear modulus of 40 Pa for this gel.

FIG. 6 :

 6 FIG.6: Scheme of a sinusoidal perturbation at the free surface of an elastic layer of height h subjected to an acceleration g. The surface at y = h is fixed.

0 e

 0 -kh a + he -kh b + e kh c + he kh d = 0 -ke -kh a + (1kh)e -kh b + ke kh c + (1 + kh)e kh d = 0.

)FIG. 7 :

 7 FIG.7: Bold line: Dimensionless acceleration (α = ρgh µ ) so that the flat interface is neutrally stable against an infinitesimal sinusoidal perturbation of dimensionless wave number k, according to Eq.16. Its minimum value α * 6.223 corresponds to the instability threshold.

2ik 2 A 2 + (s 2 +

 222 u x (h) = 0 and Eqs.3.8 are replaced by: for φ (Eq.21) and (Eq.22) in the boundaries conditions (Eqs.23-k 2 )B 2 = 0 ρg (-kA 2 + ikB 2 )ρω 2 A 1 + 2µ (k 2 A 1 -iksB 1 ) = 0 -ik (A 1 cosh kh + A 2 sinh kh)s (B 1 cosh sh + B 2 sinh sh) = 0 -k (A 1 sinh kh + A 2 cosh kh) + ik (B 1 sinh sh + B 2 cosh sh) = 0 (25)

  cosh k i k sinh k s cosh s -s sinh s k sinh k k cosh k -i k sinh s i k cosh s = 0

  k sinh k s cosh s -s sinh sk sinh k k cosh kk sinh s k cosh s kk sinh k s cos s s sin s k sinh k k cosh kk sin s k cos s = 0(27)Eqs.26 and 27 together define the dispersion relation of elastic surface waves propagating in a heavy material of height h. Plotting ω 2 as a function of k evidences that• (i) for -∞ < α < 4.5, ω 2 (k) is an increasing function, • (ii) for 4.5 • • • < α < 6.223 • • • , ω 2 (k) has a local minimum that is positive, • (iii) for α = 6.223 • • • , ω 2 (k)vanishes for a given value of k. This corresponds to the instability threshold,• (iv) for α > 6.223 • • • , ω 2 (k) is negative in a range of k. The related modes are unstable and -ω 2 (h) is their growth rate.

  1

		ρg µk -γk µ -2 e -kh	0 khe -kh	ρg µk -γk µ + 2 e kh	0 khe kh	= 0,
		-e -kh		(1 -kh)e -kh	e kh	(1 + kh)e kh
	which simplifies in
	γ µh	-e 2 k + 4 k + e -2 k	k+	ρgh µ	1 k -2 e 2 k-4

Supplementary material

A-. MEASUREMENT OF THE SHEAR MODULUS OF SOFT POLYACRYLAMIDE GELS.

An accurate measurement of the shear modulus of the materials is crucial for a quantitative analysis of the experimental observations. The linear and non linear elastic behaviours of the gels have been checked by classical rheological tests in the cone and plate geometry (using a ARES rheometer, TA Instrument). The storage modulus is much larger than the loss modulus for the relevant timescales of the experiments of the article, showing the elastic linear behaviour of this material (Fig. 5-a). Furthermore, steady stress-sweep tests evidence that the stress is proportional to the strain over a wide range of strains (Fig. 5-b), in agreement with the neo-Hookean model [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF].

In situ measurements of the shear moduli are made by indentation tests (Fig. 5-c). In order to remove surface tension contributions [3,[START_REF] Mora | [END_REF], they are done at the gel-water interface, taking care to proceed fast enough to avoid altering the mechanical properties of the gels due to water diffusion.

B-. ENERGY MINIMIZATION INCLUDING GRAVITY, ELASTICITY AND CAPILLARITY

We consider an infinite and heavy elastic layer of height h rigidly fixed under a rigid substrate.

The material is assumed to be homogeneous, incompressible, isotropic, with the density ρ and the shear modulus µ. It undergoes an outgoing acceleration g. The solid-fluid capillary constant is γ, a physical parameter that is expected to play a role if the order of magnitude of the elasto-capillary length, defined as the ratio of the capillary constant to the shear modulus, is not negligible compared to the other length-scale of the problem, i.e. the thickness h [3,[START_REF] Mora | [END_REF][37][38].

Here we show that the total energy of this system is stationary against a sinusoidal perturbation of wave number k provided the ratio ρg/µ is larger than a given value that depends on γ and k. We analytically deduce a surface tension dependant threshold for the Elastic Rayleigh Taylor Instability and the corresponding wave number.

We limit our exposition to the 2D case. A deformation is characterized by a map from the undisturbed state with coordinates (x, y) to a disturbed state R(x, y) = (X(x, y), Y (x, y)) (bold-the instability threshold.

Note that surface tension effects should be observable with thinner samples. For instance, consider a sample of thickness h ∼ 0.65 cm (made of the same gel with ρ = 10 3 kg/m 3 , g = 9.8 m.s -2 , γ ∼ 70 mN/m). Taking ρgh/µ = α * yields µ ∼ 10 and γ µh ∼ 1. Fig. 7 indicates that the threshold is shifted by the surface tension. This effect should also be observed with harder samples subjected to larger accelerations or with larger density contrasts.

C-. DISPERSION RELATION

We calculate the dispersion relation of Rayleigh waves propagating at the surface of an infinite and heavy elastic layer of height h, and rigidly fixed under a rigid substrate. The material is assumed to be homogeneous, incompressible, isotropic, with the density ρ and the shear modulus µ. It is subjected to an outgoing acceleration g. The surface tension is here neglected. We consider here plane waves with infinitely small amplitudes and the model of Hookean elasticity is used.

The propagation of elastic surface waves in heavy elastic materials has been addressed by several authors [31,32,[START_REF] Kuipers | [END_REF]40]. These authors had likely in mind the issue of waves propagating on the surface of the Earth. Hence, they were exclusively interested in accelerations directed towards the elastic medium. The effects of gravity on the wave propagation are very small in this case.

Here, the calculation of [START_REF] Kuipers | [END_REF] is modified in accordance with our experimental setup (outward gravity and no-slip condition, Fig. 6). The changes in the calculation are highlighted below.

The deformation of the elastic medium is characterized by a map from the undisturbed state (flat surface) with coordinates (x,y,z) to a disturbed state R(x, y, z) = (X(x, y, z), Y (x, y, z), Z(x, y, z)). In the following, we consider a two-dimensional problem, which amounts to imposing Z = z. The displacement u = Rr is assumed to be small, and its Cartesian components are noted u x (in plane) and u y (out of plane). Applying the decomposition


Liujiang Kang

Xiaoning Zhu

Huijun Sun

Jakob Puchinger

Mario Ruthmair

Bin Hu

Modeling the first train timetabling problem with minimal missed trains and synchronization time differences in subway networks

Keywords: First train timetabling, MIP model, Linearization, Heuristic, Case study

Urban railway transportation organization is a systematic activity that is usually composed of several stages, including network design, line planning, timetabling, rolling stock and staffing. In this paper, we study the optimization of first train timetables for an urban railway network that focuses on designing convenient and smooth timetables for morning passengers. We propose a mixed integer programming (MIP) model for minimizing train arrival time differences and the number of missed trains, i.e., the number of trains without transfers within a reasonable time at interchange stations as an alternative to minimize passenger transfer waiting times. This is interesting from the operator's point of view, and we show that both criteria are equivalent. Starting from an intuitive model for the first train transfer problem, we then linearize the non-linear constraints by utilizing problem specific knowledge.

In addition, a local search algorithm is developed to solve the timetabling problem. Through computational experiments involving the Beijing subway system, we demonstrate the computational efficiency of the exact model and the heuristic approach. Finally, three practical suggestions are proposed for the operation and management of the urban railway transit system.

Introduction

Motivation

Urban railway transit is playing an increasingly important role in reducing traffic congestion and improving the environment. One of the major roles of metro schedulers is to create timetables for the subway lines of a given network, which is a key component of railway operations and management.

According to [START_REF] Ceder | Creating bus timetables with maximal synchronization[END_REF], there are three levels of problems that have to be addressed before the actual scheduling and design of timetables can commence: (i) select the type of headway, i.e., periodic or non-periodic timetables, (ii) optimize the departure times of lines and (iii) propose objective functions, such as minimizing operator cost while providing high quality service and minimizing operator and user costs through weighting factors. Then, line operational parameters are optimized to create the whole network timetable. In addition, to comply with variations in travel demand, the day's service can be divided into several time segments (e.g., first train hours, last train hours, peak hours, etc.), where distinct timetabling strategies will be applied [START_REF] Liebchen | The first optimized railway timetable in practice[END_REF]. Table 1 presents a snapshot of the first train timetable in the Beijing subway system. The first column depicts two different transfer stations, GongZhuFen and GuoMao. The second column shows transfer directions and arrival times of the feeder trains. The third column introduces the necessary transfer times between any two lines. The fourth column provides the departure times of the connecting trains. Based on columns 1-4, the calculated transfer waiting times are listed in the final column of Table 1. The Beijing subway network example reveals that the connecting times for first trains is a serious problem. For example, if a passenger takes the first train in the up-train direction of line 1 (L1U)

Table 1 A snapshot of the first train timetable in the

and transfers to the down-train direction of line 10 (L10D) at the GongZhuFen station, a transfer waiting time of 70 min occurs. From another perspective, eight trains on L1U have departed before passengers from L10D arrive. As Figure 1 shows, first train transfers on L10D will board the 9th train on L1U. These trains are denoted as "missed trains" in this paper (Formulations are shown in Section 2).

Line departure times

Basic Data Impact factors

Type of headway

Objective function

Network Timetable

Segment running times

Station dwell times

Figure 2 Optimizing first train departure times

The total waiting time for first train transfers in the entire Beijing subway network (year 2013, 240 transfer directions) was more than 2600 min. Note that this value takes into account only the train schedule and does not consider the passengers. Therefore, this paper will study the first train timetabling synchronization problem to provide solutions that will minimize the total transfer waiting time of passengers and minimize the number of missed trains by synchronizing train arrival times at transfer stations, which will allow for well-timed connection transfers. An effective method to avoid the long-waiting-time problem as presented in this paper is depicted in Figure 2. First we collect the operational parameters (segment running times and station dwell times) according to the existing train timetable. Then, the headway and the objective function are determined. Note that under the passenger arrival pattern at stations following uniform distributions, a schedule with a constant headway between consecutive trains can reduce the total waiting time [START_REF] Niu | Optimizing urban rail timetable under time-dependent demand and oversaturated conditions[END_REF]. In addition, periodic timetables also have the advantage of being easy for passengers to memorize [START_REF] Barrena | Single-line rail rapid transit timetabling under dynamic passenger demand[END_REF].

Therefore, we adopt the fixed headway for the Beijing subway system. Finally, first train departure times in the subway network are optimized to minimize the total transfer waiting time of passengers and minimize the number of missed trains. The next section includes a review of the relevant literature before discussing the contributions of this study.

Literature review

There is a wealth of literature on periodic/non-periodic timetabling problems. [START_REF] Cordeau | A survey of optimization models for train routing and scheduling[END_REF] presented a survey of optimization models for the most commonly studied train routing and scheduling problems. [START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF] The following papers considered timetable synchronization in their studies. Ibarra-Rojas and Rios-Solis (2012) formulated a timetabling problem that maximizes the number of synchronizations in a large netwoek to facilitate passenger transfers and avoid bus bunching. the number of missed trains between line l and line l at station s .

First feeder-connecting time (FFCT) and transfer waiting time (TWT)

Before modeling the first train timetabling problem, an assumption about the supply of train services is needed. After the first train departs from the vehicle depot, we assume that the following trains in the same line are put into operation continuously with a fixed headway. Therefore, four different situations can be identified for first train transfers as shown in Figure 3. Passengers of the first feeder train (FFT) transfer either to the first connecting train (FCT), to the second connecting train or to one of the following connecting trains. In each situation, two types of transfer indicators (FFCT and TWT) are defined.

           0 D A Tra sl sl sll l t t T H         0 D A Tra l sl sl sll H t t T D A Tra sl sl sll l t t T H       
      FFCT D A Tra sll sl sl sll t t t T (1)
The TWT (which is always positive) is determined by the departure time of the connecting train and the arrival time of the first feeder passengers, see Equation ( 2), where l H  represents the headway of  l , and  sll n is an integer variable. The value of  sll n in Equation ( 3) denotes the number of trains of line  l that have passed by s when the FFT of line l arrives at station s , where function

() ceil A rounds the elements of

A to the nearest integer greater than or equal to A .

          TWT D A Tra sll sl sl sll sll l t t t T n H (2)   0, 0 | | , otherwise                D A Tra sl sl sll sll D A Tra sl sl sll l t t T n ceil t t T H (3)
Case 1:

0 TWT FFCT 
As illustrated in Figure 3 (situations 1-2), the FFT arrives at s earlier than the FCT. Therefore, passengers on the FFT board the FCT successfully. Such a case can be further divided into two situations, i.e., long waiting time (situation 1) and short waiting time (situation 2).

Case 2:

    sll l TWT FFCT n H
The FFT arrives at s later than the FCT (situations 3-4). Consequently, passengers on the FFT board the following - sll n th train, where  sll n is defined by Equation (3).

According to Figure 3 and Equations (1-3), the mathematical relationship between FFCT and TWT is presented in Equation ( 4). The first train TWT can be easily calculated as long as we know the value of FFCT. Moreover, TWT can be fixed when the arrival and departure times of the first feeder and connecting trains are obtained.

       TWT FFCT sll sll sll l t t n H (4)

First train timetable synchronization

As situation 1 in Figure 3 indicates, if the FFT arrives much earlier than the FCT then the feeder passengers have to wait a long time. Similar to situation 4, if the FFT arrives too late, the passengers on FCT must wait a long time to transfer to line l .

Time Arrival

Transfer to Transfer to To allow passengers to inter-transfer smoothly (transfer waiting time within one headway)

r a T TWT sll sll T t    Line l  Line l 1 T  1 T 2 T  2 T  l H Arrival Departure Departure ra Dw T TWT Dw l s sl l sl l ls T T t T        TWT sl l t  r a T sl l T  r a T sll T  D w l s T  TWT sll t  D w l s T ra ra TWT TWT T T Dw Dw sl l sll sl l sll l s ls l t t T T T T H             TWT TWT sl l sll l t t H      r a r a 0 0 T Dw sll l s T Dw sl l ls T T T T            
between the FFT and the FCT, the TWT is minimized:

( ) ( ) min         TWT sll sll s S l S l l L l L T v t , (6)
where  sll v represents the number of first train transfer passengers from line l to  l at station s .

Note that a different volume and spatial distribution of the first train transfer passengers in the network will lead to a different timetable and total TWT. In this paper, values of

 sll v are constant.

New first train timetabling model

As tH . From the perspective of minimum TWT, the objective function [START_REF] Chevrier | Energy saving in railway timetabling: A bi-objective evolutionary approach for computing alternative running times[END_REF] minimizes the number of trains that have passed through station s of line  l when passengers from the first train of line l arrive at s . In addition, objective [START_REF] Chevrier | Energy saving in railway timetabling: A bi-objective evolutionary approach for computing alternative running times[END_REF] ensures that the FFT and the FCT synchronize at transfer stations.

( ) ( ) ( ) ( ) min ( ) ( ) ( )                                   AA sll sll l sl l sl l l sll sl l sl sl s S l S l l L l L s S l S l l L l L N v n H v n H v v t t (7)
Remark 1: The total TWT of optimal solutions obtained by objective function [START_REF] Chevrier | Energy saving in railway timetabling: A bi-objective evolutionary approach for computing alternative running times[END_REF] for first train transfers is equivalent to that of objective function [START_REF] Ceder | Creating bus timetables with maximal synchronization[END_REF].

Proof: For any lL  , lL   , ( ) ( )   s S l S
l , and inter-transfer directions (   ll and

  ll ), we have           TWT D A Tra sll sl sl sll sll l t t t T n H and          TWT D A Tra sl l sl sl sl l sl l l t t t T n H . Therefore,        TWT TWT sll sll sl l sl l v t v t (8) ( ) ( )
                     D A Tra D A Tra sll sl sl sll sll l sl l sl sl sl l sl l l v t t T n H v t t T n H ( ) ( )                         A Dw A Tra A Dw A Tra sll sl l s sl sll sll l sl l sl ls sl sl l sl l l v t T t T n H v t T t T n H ( ) ( )                                  A A
v n H v n H v v t t v T v T v T v T ( ) ( )                    AA sll sll l sl l sl l l sll sl l sl sl v n H v n H v v t t C ,
where passenger volumes

 sll v and  sl l v are given based on real data. Polynomial                Dw Dw Tra Tra sll l s sl l ls sll sll sl l sl l C v T v T v T v T is a constant. Note that the headways ( l H
) for the first trains in a subway network are given based on different transit lines. Therefore, the total first train TWT in a subway network can be captured by Equation ( 9), ( ) ( )

         TWT sll sll s S l S l l L l L vt (9) ( ) ( ) 1 () 2                TWT TWT sll sll sl l sl l s S l S l l L l L v t v t   ( ) ( ) 1 ( ) ( ) 2                           AA sll sll l sl l sl l l sll sl l sl sl s S l S l l L l L v n H v n H v v t t C ( ) ( ) ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) 2                                                  AA sll sll l sl l sl l l sll sl l sl sl s S l S l l L l L s S l S l l L l L s S l S l l L l L C v n H v n H v v t t ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) 2                                        AA sll sll l sl l sl l l sll sl l sl sl s S l S l l L l L s S l S l l L l L v n H v n H v v t t , where ( ) ( ) 1 2         s S l S l l L l L C is a constant.
Therefore, objective function [START_REF] Chevrier | Energy saving in railway timetabling: A bi-objective evolutionary approach for computing alternative running times[END_REF]

T T T (13)

According to the real data from the Beijing subway system, dwell time in this paper lasts 30 s at an ordinary station and reaches 60 s at a transfer station.

Line headways in the Beijing subway depend on the time period. During rush hours (7:00 to 9:00 and 17:00 to 19:00), the headway is 2 min. Generally, the headway during non-rush hour periods is 5 min. The headway for the first/last trains is set to a maximum of 10 min.

First-train timetabling models

The model M1 in Equation ( 14) is a non-linear model and aims to minimize the total TWT for first train transfers (with passenger volumes

: min M1         TWT sll sll s S l S l l L l L T v t (14) . . , , , ( ), ( )
                 TWT D A Tra sll sl sl sll sll l s t t t t T n H l L l L s S l s S l 11 , , (
)          ss D R Dw sl l lk lk kk t t T T l L s S l 11 , , ( )
              ss D R Dw sl l l k l k kk t t T T l L s S l 1 11 , , ( )
          ss A R Dw sl l lk lk kk t t T T l L s S l 1 11 , , ( ) 
              ss A R Dw sl l l k l k kk t t T T l L s S l   0, 0 , , , ( ), ( ) |
| , otherwise                       D A Tra sl sl sll sll D A Tra sl sl sll l t t T n l L l L s S l s S l ceil t t T H int tra , { , }, , , ( ), ( )
         Dw Dw Dw Dw ls l s T T T T l L l L s S l s S l min max ,     l l l T t T l L min max ,         l l l T t T l L
There are non-linear constraints (for determining  

0, 0 0 | | , otherwise                          D A
: min ( ) ( ) M2                                 AA sll sll sl l sl l sll sl sl s S l S l l L l L s S l S l l L l L N n n t t (17) . . , , , ( ), ( )
             sll sll l s t v H l L l L s S l s S l , , , ( ), ( )
           sl l sl l l v H l L l L s S l s S l , , , ( ), ( )
             sll sll sl l v v l L l L s S l s S l 11 , , (
)          ss D R Dw sl l lk lk kk t t T T l L s S l 11 , , ( )
              ss D R Dw sl l l k l k kk t t T T l L s S l 1 11 , , ( )
          ss A R Dw sl l lk lk kk t t T T l L s S l 1 11 , , ( )
               ss A R Dw sl l l k l k kk t t T T l L s S l 0, , , ( ), ( )
                D A Tra sl sl sll sll l t t T n H l L l L s S l s S l int tra , { , }, , , ( ), ( )
         Dw Dw Dw Dw ls l s T T T T l L l L s S l s S l min max ,     l l l T t T l L min max ,         l l l T t T l L
3 Solution algorithm and sample tests

Local search heuristic

The number of integer variables in the above integer programming model is a good index of its complexity. Variable  sll n represents the missed trains between the first train of line l and the first train of line  l at station s . This means that an integer variable is used for every combination of two first trains on different lines that intersect at s . The number of integer variables in the worst case is that would solve such problems in a reasonable time due to its scalability.

l t 1 l t  1 l t  Original solution l l t t     New solution 1 l t  l t 1 l t  Solution expression   Line Number l 1 l  1 l    Time l t  1 l t  1 l t  Time Line Line

Figure 5 Producing new solutions

The heuristic (Algorithm 1) contains an outer loop to select a sensitive metro line and an inner loop to optimize its departure time. We use d is copied to d  , and operators in Figure 5 are applied to d  to produce new solutions. In Figure 5, each solution is coded in a vector according to line numbers. We select a number l according to Algorithm 2 and take solution operator to l t (one-to-one correspondence). The original departure time l t of line l is replaced by

  l t
, where parameter  is randomly generated as a mutation factor and is limited with Equation [START_REF] Ibarra-Rojas | Synchronization of bus timetabling[END_REF]. A new solution is created after updating the elements in the vector. Then, we use the methods shown in Figure 6 to make timetables for solutions best d and d  . In Figure 6, based on existing running times and dwell times, the (network) first train timetable can be produced with constraints [START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF][START_REF] Ibarra-Rojas | An integrated approach for timetabling and vehicle scheduling problems to analyze the trade-off between level of service and operating coasts of transit networks[END_REF][START_REF] Kang | A case study on the coordination of last trains for the Beijing subway network[END_REF] when the line departure times are calculated from the above step.

Further, the periodic timetable can also be designed with the proper headway included. According to the above timetables, the corresponding objective functions are calculated. Following that, the algorithm decides whether to accept the new solution or not. After that, the inner loop p I and the out loop i O are updated. Finally, the algorithm stops when the best solution has been found.

Algorithm 1.

Local search heuristic

1:

Initialize: read running times and dwell times of the original timetable; input initial solutions.

2:

best d d initial d  .

3:

For max 1  i OO 4:
Search for the most sensitive line l using Algorithm 2.

5:

For

7:

Generate new solution from the above step, dd   .

8:

Create timetables for d  and d // See Figure 6.

9:

If ( ) ( )   N d N d ,

2:

Change line departure time ii tt , ii tt .

3:

Update the upper and lower departure times in mm  matrixes D and D .

4:

End for

5:

With the solutions in matrixes D and D , create timetables according to model [START_REF] Liebchen | The first optimized railway timetable in practice[END_REF].

6:

Calculate the objective values for the above timetables.

7:

Store the objective values in two 1 m matrixes ( () ND and (D) N ), Equations [START_REF] Liebchen | Computing delay resistant railway timetables[END_REF][START_REF] Lin | Using genetic algorithms to optimize stopping patterns for passenger rail transportation[END_REF].

8:

Select the largest value by ( ) (D)

1 1 1 1 1 1 1                      i i i i i i i m m m m m t t
(18) 1 1 1 1 1 1 1                              i i i i i i i m m m m m t t
| ( ) (D) |  ll N D N ) in   ,   ll tt
. However, parameter  should be limited to as small a number as possible under the condition of reasonable algorithm efficiency. □

Sample tests

A sample network with three bi-directional lines crossing at two transfer stations is designed for the numerical tests; see Figure 8 for an illustration. Segment running times between two stations are given along the lines. In addition, the dwell time at any transfer station is set to 1 min ( ). However, 12 connections force passengers to wait for more than one headway time. Among them, four connections in Table 3 have as large as three missed trains, which can also be improved. As for the TWT, passengers wait 1605 min in total with the original timetable.

5 R l s T  1 0 R ls T  5 R l s T  5 R l s T  4 R l s T  4 R l s T  4 R l s T  1
-- -- -- -- -- -- 20 1605
The optimal results in Table 4 are found by solving the MIP model, using the heuristic algorithm and the CPLEX. Detailed analyses of optimization results are provided below.

The number of missed trains is reduced from 20 to 8, an improvement of 60.0%. A common phenomenon in the "  sll n " column of Table 4 is that there is only one inter-transfer direction at most being synchronized by two first trains. To explain this phenomenon, we take line 1U and line 2U crossing at station A in Figure 8 as an example. If both directions are connected by the first trains at station A , then we will have Equation (21) and Equation (22). However, Equation ( 22) is not true because we set 3   TT in this test. Therefore, there is only one direction of two inter-transfer directions that can be synchronized by two first trains.

1U2U 2U 2 U 1U2U 2 2U 1U2U 2U 2U 1U 1U 1U 1U 2U 2U 2 1U 1U 1U 2U1U 1 U 2 U U 1 U1 0 0                               Tra A Dw Tra Dw A Tra A A A A A A A Tra A Dw Tra Dw A Tra A A A A A A A D A A A A A A A D A A A A A A A T t T T T t T T t T t t t t t t t T T T t t (21) 1U 2U 2U1U 1U2U    Dw Dw Tra Tra A A A A T T T T ( 22
)
The total TWT has been reduced from 1605 min to 345 min, a decrease of 78.5%. As Table 4 shows, passengers can transfer to the FCT ( 0   sll n

) or the second train (

1   sll n ) of each line
smoothly. In addition, the TWT of each direction is controlled within reasonable limits.

In addition, long waiting times of more than one headway have been completely eliminated. If the first feeder train arrives earlier than the first connecting train for one headway time, then the first feeder passengers will have a long wait. Therefore, the optimal timetable coordinates the arrival/departure times of the first feeder and the connecting trains to reduce the TWT for passengers effectively. where  is limited to [-5, 5] min, the six first trains depart from their depots at 4:56:00, 4:57:00, 5:04:00, 5:05:00, 5:04:00 and 5:05:00. The performance of the heuristic algorithm is compared with that of CPLEX. In this small example, the heuristic algorithm takes less than 1 s and is always able to find the optimal solution, whereas CPLEX takes 18 s to return a solution. As illustrated, the arrival times of the first feeder trains, the transfer times and the departure times of the first connecting trains are obtained from the published first train timetable. The number of missed trains and the corresponding TWT are calculated and listed in the last two columns of Table A.

Optimal Beijing subway timetable

The MIP model was solved using CPLEX 12.6.1 on an Intel Core i5 PC with 2.5 GHz and 6 GB RAM. The optimal first train timetables for the Beijing subway system, which are optimized under departure time constraints ranging from 4:50:00 to 5:10:00, reduce the number of missed trains from 441 to 424 and the total TWT from 22790 min to 16380 min. The optimization results of the MIP model with respect to line 1 are presented in Table B (see Appendix). In this case, adjusting the arrival and departure times of the first feeder and connecting trains reduced the number of missed trains from 85 to 79. This results in a total TWT of 6774 min, a decrease of 19.8% from the original value.

Figure 11 shows the time-distance diagrams of each line at several representative transfer stations in this network. It is unavoidable that some first train transfer passengers will have long wait times in a large subway network due to the network structure, such as at the GongZhuFen station on line 1 and line 10. An effective method to diminish this problem is to stagger line departure times. Another solution is for subway companies to operate half-way trains to eliminate departure time differences. Figure 12 shows the distributions of waiting times for all transfer directions. As illustrated, the optimal timetables increase the number of directions which have less TWT, e.g., TWT in [0,1] min and (1,5] min. With respect to long TWT (TWT in (5,10] min, (10,30] min and (30,70] min), the optimal timetables reduce the number of directions. In our proposed first train model, departure times are limited between 4:50:00 and 5:10:00 to ensure that all lines begin operations around the same time. However, it is interesting to test the MIP model without use of Equation [START_REF] Ibarra-Rojas | Synchronization of bus timetabling[END_REF], which sets constraints on line departure times. Without departure time limits, first trains can start service at any time of the day (theoretically from 00:00 to 24:00). The results in Table 5 compare the above two cases. First, we use the CPLEX to solve the MIP model with first train departure time limits (with Equation 12). As a result, the number of missed trains in the Beijing subway network is 424 and the TWT reaches 16380 min. When solving the model without Equation ( 12), the number of missed trains is 408 and the TWT is 13550 min. We then use the heuristic to solve the model again. The results show that the number of missed trains is 429 and the TWT is 16690 min when using Equation [START_REF] Ibarra-Rojas | Synchronization of bus timetabling[END_REF]. Meanwhile, the number of missed trains reduces to 419 and the TWT descends to 13972 min when Equation ( 12) is not applied. In addition, the heuristic only takes approximately 12 s to run. Therefore, the heuristic has potential to solve even larger networks, such as the London Network, the New York Network and the Tokyo Network.

Regarding the effect of an initial solution, we also tested the heuristic with and without an initial solution. No significant differences in the solutions returned were observed. Therefore, this heuristic

Figure 13 Comparisons of solutions obtained with/without departure time limits

To compare the above two optimal solutions intuitively, a bar chart of different decision variables is provided in Figure 13, where the x-axis shows 18 operating lines and the y-axis represents the first train departure times. As illustrated, most lines have similar departure times with or without the start time limits, with the exception of lines 9D and 10U. Line 9D departs at 5:16:00 with limits and departs at 5:37:00 without limits. Line 10U departs at 4:55:00 with limits and departs at 4:35:00 without limits.

These results are also acceptable in real operations. Compared to the original departure times, most lines have similar decisions, with the exception of line 9D, which is due to particular factors that were considered in real operations (low passenger flows in the early morning in line 9D).

Effects of train dwell times on TWT

As Equations [START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF][START_REF] Ibarra-Rojas | An integrated approach for timetabling and vehicle scheduling problems to analyze the trade-off between level of service and operating coasts of transit networks[END_REF] show, the arrival and departure times are influenced by dwell times. In addition, Equations indicate that train dwell times play an important role in timetable synchronization. Therefore, we test the effects of train dwell times on the solutions in Table 6. The case of the 60-s dwell time at transfer stations was selected as the benchmark and results in 424 missed trains with 16380 min of total TWT. The number of transfer directions of TWT within [0,1] min, [START_REF] Barrena | Single-line rail rapid transit timetabling under dynamic passenger demand[END_REF][START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF] min, [START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF][START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF] min and (10, ]

 min are 58, 75, 52 and 55, respectively. This experiment indicates that large dwell times can reduce missed trains, but they also increase TWT for passengers.

L1U L1D L2U L2D L4U L4D L5U L5D L9U L9D L10U L10D L6U L6D L8U L8D L13U

Conclusions

The problem of first train timetabling considering passenger transfers in a large subway network is an inevitable and challenging optimization problem. Planning departure times of subway lines reasonably can decrease missed trains as well as passenger waiting times at transfer stations. This paper proposes a new approach to precisely describe the first train transfer problem in mathematical terms.

Then, non-linear constraints are linearized. Through suitable adjustments of arrival/departure times of feeder/connecting trains, the number of missed trains and passenger waiting times for first train transfers can be minimized. Computational experiments have been carried out on the Beijing subway system based on the original timetables.

Based on research results, there are three general suggestions put forward for the operations and management of urban railway transit. First, headways of the first trains can be reasonably increased to the maximum values. This not only saves operating costs for companies but also reduces the number of missed trains. From the passengers' perspective, transfer waiting times do not increase. In addition, operating strategies such as running interval trains can be applied to large loop lines and long distance lines to decrease waiting times for first train transfers. Taking Beijing subway line 10 for instance, the whole trip time is 133 min. With interval trains, the long waiting times listed in Table 1 can be avoided effectively. Finally, the use of two or more vehicle depots could significantly help to reduce transfer waiting times. They can be located on the diagonals of loop lines or in the middle of straight lines.

The model proposed can be expanded in several ways. For example, lines, transfer stations and directions can be defined with weights, which would make it a weighted optimization problem. In this paper, we assign passengers to different transfer directions to make the model more practical. In addition, the first train timetabling model could be extended to a more comprehensive period-based timetabling model for all the transition periods, from the first train to the morning peak hours. In this way, train intervals of each line can be reduced gradually through inserting additional trains according to the increase in transportation demands. Other methods such as stochastic theory and real-time dispatching can also be considered for further study.

Figure 1

 1 Figure 1 Affiliated problem of long transfer waiting time

  presented a global review of the crucial strategic and tactical steps of the design and scheduling of the network. Cacchiani et al. (2014) presented an overview of real-time railway rescheduling. When designing public transport timetables, it is common to optimize an objective function relevant for the operating companies, such as minimizing the operating costs (Ibarra-Rojas et al. 2014), minimizing the train/vehicle running times (Cacchiani and Toth, 2012), maximizing the utilization of energy (Li and Lo, 2014) or maximizing the profit earned by the company (Lin and Ku, 2014). Chevrier et al. (2013) considered the optimization of the running times, a step in train planning that precedes the timetabling. The calculation was directly related to the construction of speed profiles, which indicate the speed that the trains must hold at each position. Using a directed multi-graph, Caprara et al. (2002) proposed a graph-theoretic formulation for the train timetabling problem in which nodes correspond to departures/arrivals at a certain station at a given time. The formulation was used to derive an integer linear programming model that used the Lagrangian relaxation method. Niu and Zhou (2013) constructed a binary integer programming model incorporated with passenger loading and departure events to optimize a passenger train timetable in a heavily congested urban rail corridor. Fischetti et al. (2009) proposed and computationally analyzed four different methods to improve the robustness of a given train timetabling problem solution for the non-periodic case. Liebchen et al. (2010) studied the design of delay-resistant periodic timetables to optimize the transfer time between two adjacent trains to not only ensure successful transfers but also to minimize total travelling cost. Real-time railway disruption management is currently an active area of study in Operations Research, and it includes real-time timetable rescheduling and real-time rescheduling of the rolling stock (Cacchiani et al., 2014). From the infrastructure point of view, one widely proposed objective is to minimize the deviation from a solution proposed by the operator. Kang et al. (2015b) minimized the difference between the published timetable and the rescheduled one with train delay conditions. Pellegrini et al. (2014) proposed a mixed-integer linear programming formulation for tackling train routing and scheduling in case of disruption. Yin et al. (2016) proposed a stochastic programming model for the metro train rescheduling problem, which was reformulated into a decision-making process. The goal of Cacchiani et al. (2016) was to design conflict-free schedules that differ as little as possible from the ideal ones. They proposed an integer linear programming model for this problem, and the linear programming relaxation of the model was used to derive a dual bound. Because of disruptions, real-time scheduling and routing decisions are frequently made. Lamorgese et al. (2016) presented a mixed integer linear programming model to minimize deviations from the official timetable, and solved it with a Benders'-like decomposition within a suitable master scheme. In a passenger railroad system, Lin and Ku (2014) addressed the stopping pattern optimization problem to determine the train stopping strategy to maximize the profit made by a rail company. Liu et al. (2013) proposed three objectives to reflect the benefits of bus operator and/or passengers, including minimizing the total waiting time, the total in-vehicle travel time and the total operating cost. Li and Lo (2014) formulated an integrated energy-efficient operation model to jointly optimize the timetable and speed profile to achieve a better performance on the net energy consumption (the difference between the tractive energy consumption and the utilization of regenerative energy). Considering the uncertain train mass, variable tractive force, braking force, and running resistance, Yang et al. (2016) developed a stochastic programming model to determine the timetable and speed profile for the subway system. A simulation-based genetic algorithm incorporated with the optimal control algorithm was designed to solve the model. Several papers optimize an objective function relevant to passenger waiting time. Kang et al. (2015a) used a mathematical method to reveal the relationships between passenger transfer connection time and passenger transfer waiting time of last subway trains. Wong et al. (2008) presented a mixed-integer programming optimization model for the train synchronization problem for non-periodic timetables that minimized the transfer waiting times of all passengers. Shafhi and Khani (2010) formulated a mixed-integer programming model that gave the departure times of vehicles in lines so that passengers could transfer between lines at transfer stations with minimum waiting times. Niu et al. (2015) studied the train scheduling problem that considers time-dependent demand and skip-stop pattern, and proposed non-linear integer programming models to minimize passenger waiting times.

Figure 3

 3 Figure 3 Relationships of FFCT and TWTThe FFCT is influenced by the departure time of the first connecting train and the arrival time of the first feeder passengers, see Equation[START_REF] Barrena | Single-line rail rapid transit timetabling under dynamic passenger demand[END_REF], where

Figure 4 H

 4 Figure 4 First train timetable synchronization in an ideal situation Figure 4 shows how two first trains in different lines synchronize at a transfer station in an ideal situation. The time difference between 1T departure and 1 T arrival is

  mentioned above, two indicators, FFCT and TWT, are proposed to evaluate the first train transfers. If the FFT and the FCT synchronize at s successfully, then 0 Otherwise, one side of transfers may face a long waiting time, i.e.,



  sll v representing the number of first train transfer passengers from line l to  l at station s ) in a subway network. Model M1 was adopted from Wong et al. (2008) and Shafhi and Khani (2010) to calculate the total passenger transfer waiting time in a given network. ( ) ( )



  sll n ) in model M1 that make the model difficult to solve by mixed integer programming (MIP) solvers. Therefore,  lL , are reformulated as follows.

S

  denotes the number of transfer stations and L denotes the number of lines in the network. If we want to solve a large network, e.g., the Pairs network or the London metro, it would require much time. This was the motivation for the development of an effective local search algorithm

O

  to record the number of loops that run in the inner and outer loop algorithms, respectively. First, an initial solution d (randomly produced or input) is set as the best solution best d . The values of max I , max O ,

Figure 6 2 :m

 62 Figure 6 Creating timetables for the first trainsAlgorithm 2: Sensitivity test of line departure times in the network. This method selects one of the most sensitive lines (the largest absolute value difference of objective functions between the upper and lower bounds of associated solutions) to search for the local best-case scenario. As Figure7shows, the

Figure 7

 7 Figure 7 Modifying line departure times

Figure 8 A

 8 Figure 8 A sample network for numerical tests

Figure 9

 9 Figure 9 Optimal line departure times of the sample network Departure time variations of lines can be found in Figure 9. As mentioned before, the original line departure times (

Figure 10

 10 Figure 10 Map of the Beijing SubwayTable A in the appendix shows the results of missed trains and transfer waiting times of 56 transfer directions for the original Beijing subway system. For space reasons, we only show the data for Line 1.

Figure 11

 11 Figure 11 Time-distance diagrams at representative transfer stations

Figure 12 TWT

 12 Figure 12 TWT bar chart of two first train timetables for the Beijing subway 4.3 With and without constraints of first train departure times

  1] min[START_REF] Barrena | Single-line rail rapid transit timetabling under dynamic passenger demand[END_REF][START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF] min[START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF][START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF] min[START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF]30] min (30,70] minDirectionsOriginal Optimal algorithm can be implemented to solve other timetabling problems and is quite robust because its performance is not influenced by the quality of the initial solution[START_REF] Barrena | Single-line rail rapid transit timetabling under dynamic passenger demand[END_REF]).

Table 2 . Table 2 Notations used in this paper Parameters and constants:

 2 [START_REF] Ceder | Creating bus timetables with maximal synchronization[END_REF] formulated as a mixed-integer linear programming problem to address the problem of generating a timetable for a given network that maximized the number of simultaneous bus arrivals at the transfer nodes. Wu et al.

	(2015) proposed a timetable synchronization optimization model to optimize passengers' waiting time
	while limiting the waiting time equitably over all transfer stations, which improves the worst transfer in
	the subway network.
	In a previous paper, Kang et al. (2015a) studied the last train transfer problem, which focused on
	network connections. Transfer results were mathematically represented by the network accessibility
	and transfer waiting times. The main difference in approach between the first train timetabling problem
	and the last train timetabling problem is as follows. In the first train problem, transfers in a complex
	subway network often face long connecting times that have not been well addressed. This issue has not
	received enough attention by researchers yet, although it has been raised by the Beijing traffic control
	center and urban railway transit industry. However the last train problem addresses the transfer
	possibility, which is a binary problem. Therefore the mathematical models differ substantially. In this
	paper, we consider the first train transfer problem in a large subway network. Specifically, we first
	provide a new formulation to precisely describe this problem in mathematical terms. To improve the
	solution performance of the proposed model, the non-linear constraints are linearized. In addition, a
	heuristic algorithm is developed to solve the Beijing subway case with a local search method. Solutions
	obtained by this algorithm for variables, missed trains, passenger waiting time and CPU time, are
	compared with those by CPLEX to demonstrate the computational efficiency of the exact model and
	the heuristic approach.
	The remainder of this paper is organized as follows. In Section 2 we formally describe the
	first-train timetabling problem and present two first train timetabling models. In Section 3 we describe
	our algorithm and test its efficiency with a sample case. Then we present the results of the Beijing
	subway case in Section 4, and conclusions in Section 5.
	2 Model formulations
	In this section, two indicators, first feeder-connecting time (FFCT) and transfer waiting time
	(TWT), are introduced to analyze the first train transfer problem. The mathematical relationship
	between FFCT and TWT is given, and a new first train timetabling model is proposed to minimize the
	total passenger transfer waiting time and minimize the number of missed trains	 n that have passed sll

  can achieve the minimum TWT for first train transfers in subway

	systems. □
	2.4 Constraints of the first train timetabling
		This section gives the following six constraints for the first train timetabling, including arrival
	times at stations, departure times from stations, departure times from vehicle depots, segment running
	times, station dwell times and line headways.
		The first train arrival time at station s	in line
									ss	1 
				t	A sl	  l t	11 kk   R lk  T	Dw lk T	(10)
									ss
		t	D sl	  A sl t	Dw ls T	  l t	11 kk   R lk  T	Dw lk T	(11)
		The first train departure time from the vehicle depot is limited by Equation (12), where min T
									l
	represents the earliest first train service time and max T	represents the latest service time.
									l
							T	min	 t T	max	(12)
								l	l	l
		Segment running times are obtained based on the original Beijing subway timetable, as published
	in 2014.
		Stations are divided into intermediate stations and transfer stations. In general, trains dwell at
	stations for different amounts of time depending on passenger flows and various station operations,
	such as changing the arrival-departure tracks or changing train directions. In subway systems, train
	dwell times are mainly predetermined according to station types, i.e., intermediate station dwell time
	(	Dw T ) and transfer station dwell time (	Dw T ). Therefore, Dw T	is limited to the following Equation
		int					tra		ls
	(13).
							Dw		{	Dw	,	Dw	}
							ls		int	tra

l can be tracked by Equation

[START_REF] Guihaire | Transit network design and scheduling: A global review[END_REF]

, where l t denotes the departure time from the vehicle depot in line l , R ls T denotes the running time between the previous station and station s in line l , and Dw ls T denotes the dwell time at station s in line l .

Similarly, the first train departure time at station s in line l is captured by Equation

[START_REF] Ibarra-Rojas | An integrated approach for timetabling and vehicle scheduling problems to analyze the trade-off between level of service and operating coasts of transit networks[END_REF]

.

  have passed through station s in line  l when passengers arrive at s . The network based TWT for passengers can be minimized by the objectives in model M1 (or model M2). Clearly, values of TWT are greater than or equal to zero. Therefore, we have the following Equation (16), which also determines As proved previously, to reach the minimum TWT goal while taking passenger volumes into account, model M1 can be reformulated as model M2. The first train timetabling model M2 in Equation (17) aims at minimizing the number of missed trains while considering transfer passengers. On the one hand, M2 can be indirectly understood to minimize the total first train TWT for passengers in a subway network. On the other hand, this model also reduces the number of operating trains, which is important from the viewpoint of railway companies.

	Remark 2: Constraint	t	 D sl	 Tra sll t T   A sl		 n H  sll	l			0	can directly capture the correct number of
	sll missed trains that  n . With the minimization objective of M2, the
	non-linear constraints can be represented by Equation (16). □
						t	 TWT sll	    sl sl D A t t	 Tra sll T		sll n			H	l			0	(16)
	( )	( )												( )	( )
						Tra
	t		t T
	sll n	sl	sl		sll									t	D sl	A sl t T sll Tra	sll n H	l	(15)
			D	A		Tra
	ceil t		t T	H
			sl	sl		sll	l

Proof:

Model M2 represents our final MIP model and can be solved directly by any MIP software (e.g., CPLEX) or standard optimization package (e.g., GAMS).

  1  passengers is given in the second column of Table3. For comparison, each first train starts at 5:00:00 am, and the transfer results of the original timetable are recorded in Table3. As we can see, there are 20 missed trains in this sample test. First-to-first transfers are only possible for

																T	Dw	);
																ls
	transfer time between two lines is set to 3 min (	sll Tra   T	3	); the headway time of line 1 is 10 min
	( 1U HH 1D 	10	); and headways of line 2 and line 3 ( 2U H		H	2D		H	3U		H	3D		5	) are 5 min. The
	number of transfer four directions (	 	0

sll n

Table 3 First train connections of the original timetable

 3

	Direction	sll v		Feeder Line Vehicle Dep. Arrival	Transfer	Connecting Line Vehicle Dep. Departure	n sll		TWT (min)
	1U to 2U	10	5:00:00	5:05:00	3 min	5:00:00	5:06:00	1	30
	2U to 1U	30	5:00:00	5:05:00	3 min	5:00:00	5:06:00	1	240
	1U to 2D	10	5:00:00	5:05:00	3 min	5:00:00	5:05:00	1	20
	2D to 1U	10	5:00:00	5:04:00	3 min	5:00:00	5:06:00	1	90

Table 4 First train connections of the optimal timetable

 4

	Direction	Feeder Line Vehicle Dep. Arrival	Transfer	Connecting Line Vehicle Dep. Departure	sll n		TWT (min)
	1U to 2U	4:56:00	5:01:00	3 min	5:05:00	5:10:00	0		60
	2U to 1U	5:04:00	5:09:00	3 min	4:57:00	5:02:00	1		0
	1U to 2D	4:56:00	5:01:00	3 min	5:06:00	5:10:00	0		60
	2D to 1U	5:05:00	5:09:00	3 min	4:57:00	5:02:00	1		0
	1D to 2U	4:57:00	5:12:00	3 min	5:05:00	5:10:00	1		0
	2U to 1D	5:04:00	5:09:00	3 min	4:58:00	5:13:00	0		20
	1D to 2D	4:57:00	5:12:00	3 min	5:06:00	5:10:00	1		0
	2D to 1D	5:05:00	5:09:00	3 min	4:58:00	5:13:00	0		10
	1U to 3U	4:56:00	5:12:00	3 min	5:05:00	5:10:00	1		0
	3U to 1U	5:04:00	5:09:00	3 min	4:57:00	5:13:00	0		15
	1U to 3D	4:56:00	5:12:00	3 min	5:06:00	5:10:00	1		0
	3D to 1U	5:05:00	5:09:00	3 min	4:57:00	5:13:00	0		30
	1D to 3U	4:57:00	5:01:00	3 min	5:05:00	5:10:00	0		60
	3U to 1D	5:04:00	5:09:00	3 min	4:58:00	5:02:00	1		0
	1D to 3D	4:57:00	5:01:00	3 min	5:06:00	5:10:00	0		90
	3D to 1D	5:05:00	5:09:00	3 min	4:58:00	5:02:00	1		0
	Total	--	--	--	--	--	8		345

Table 5 With and without first train departure time limits

 5

	Solution algorithm	Departure time limits	Missed trains	TWT	CPU Time
	CPLEX	With Equation (12) Without Equation (12)	424 408	16380 min 13550 min	12 min 35 s 58 min 12 s
	Heuristic	With Equation (12) Without Equation (12)	429 419	16690 min 13972 min	12 s 14 s

Table 6 Analysis of train dwell times

 6

	Dwell time	,,  sll  s l l n		,,  TWT   sll s l l t	[0,1]min (1,5]min	(5,10]min (10,30]min	(30, 70]min
	30 s	424		16380	66	79	40	34	21
	60 s	424		16380	58	75	52	35	20
	90 s	422		17150	24	80	76	40	20
	120 s	318		17340	33	73	72	41	21

Acknowledgements

This paper is supported by the China National Funds for Distinguished Young Scientists (No. 71525002), the NSFC (No.71473259, No.71322102). The authors thank the editor (Prof. Hai Yang) and the referees for their meticulous review and their valuable suggestions, which greatly improved this paper. In addition, this work has been partially financed by the Austrian Federal Ministry of Science, Research and Economy (BMWFW) within the EURASIA PACIFIC UNINET scholarship framework.

Appendix Table A. Original first-train connections of Line 1 of the Beijing Subway

Station Direction

 sll v

First Feeder