Quan Jin

Claude Verdier

Pushpendra Singh

Nadine Aubry

Roxana Chotard-Ghodsnia

Alain Duperray

Nadine Aubry ++

Migration and deformation of leukocytes in pressure driven flows

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Although the white blood cells (WBC or leukocytes) constitute only a small fraction of blood cells, they play a very important role in the body's immune response. They do so either by fighting pathogens, viruses, or other cells (cancer cells for example), or by reaching the site of infection through the transendothelial migration. The endothelial monolayer is a part of the vessel wall. Endothelial cells are maintained in close contact with each other by a complex network of transmembrane adhesion proteins, especially cadherins. Diapedesis is the way leukocytes migrate through the endothelial junctions under the combined effects of signalisation, changes in the rheological properties and the cell-cell adhesion properties. As part of the migration process, circulating leukocytes must first adhere to the luminal surface of the endothelium (see figure 1).

The goal of this paper focuses on the first phase of the above process, i.e., the study of the motion and deformation of leukocytes in small vessels and the investigation of conditions under which leukocytes can remain attached to the vessel wall and do not drift away. For this, we will assume that the leukocyte is already close to the endothelium. The factors, that cause leukocytes to get close to the endothelium in the first place, such as collisions with other cells, will not be studied. This work should also be helpful in improving our understanding of the motion of cancerous cells, and the role of the flow rate and the cortical tension in the process [START_REF] Burdick | Colon carcinoma cell glycolipids, integrins, and other glycoproteins mediate adhesion to HUVECs under flow[END_REF]. Indeed, it has been observed that the mechanisms by which cancer cells [START_REF] Burdick | Colon carcinoma cell glycolipids, integrins, and other glycoproteins mediate adhesion to HUVECs under flow[END_REF][START_REF] Tözeren | E-selectin-mediated dynamic of breastand colon cancer-cells with endothelial[END_REF][START_REF] Haier | Tumor cell adhesion under hydrodynamic conditions of fluid flow[END_REF] migrate to and through the endothelium to reach different tissues are similar to the ones used by leukocytes, although the adhesion molecules involved for the former are not all known yet [START_REF] Chotard-Ghodsnia | Morphological analysis of tumor cell/endothelial cell interactions under shear flow[END_REF].

Past recent numerical studies [START_REF] Kan | Hydrodynamics of a compound drop with application to leukocyte modeling[END_REF][START_REF] N'dri | Computational modeling of cell adhesion and movement using a continuum-kinetics approach[END_REF][START_REF] Khismatullin | Three-dimensional numerical simulation of a receptor-mediated leukocyte adhesion to surface: Effects of cell deformability and viscoelasticity[END_REF], based on a compound-drop model, have been capable of qualitatively explaining the behavior of leukocytes observed in experiments. The compound-drop model of a leukocyte consists of a spherical nucleus containing a core fluid which is surrounded by a thick layer of cytoplasm. In this paper, we use both this compound-drop model and a drop-rigid-particle model to represent leukocytes. The two most common adhesion models are based on Bell's equilibrium model and Dembo's kinetics model [START_REF] Bell | Cell adhesion. Competition between nonspecific repulsion and specific bonding[END_REF][START_REF] Dembo | The reaction-limited kinetics of membrane-to-surface adhesion and detachment[END_REF][START_REF] Lauffenburger | Receptors: Models for binding, trafficking, and signaling[END_REF][START_REF] Sukumaran | Influence of shear flow on vesicles near a wall[END_REF], which are both used here to explore the adhesion mechanisms.

Figure 1. The different steps corresponding to the leukocyte extravasation.

MODELING AND GOVERNING EQUATIONS 2.1 Rheological models for leukocytes

In this study leukocytes, which are the only blood cells with a nucleus, will be represented using two different two-layered models in which the outer layer represents the cytoplasm and the inner layer represents the nucleus [START_REF] Kan | Hydrodynamics of a compound drop with application to leukocyte modeling[END_REF][START_REF] Khismatullin | Three-dimensional numerical simulation of a receptor-mediated leukocyte adhesion to surface: Effects of cell deformability and viscoelasticity[END_REF]. The two-layered models of leukocytes used here approximately represent neutrophils, although for simplicity the shape of the nucleus is assumed to be spherical.

In the first two-layered model shown in Figure 2a, called the drop-rigid-particle model, the cytoplasm is assumed to be a Newtonian fluid, and the inner core is assumed to be rigid and thus modeled as a rigid particle. The nucleus volume is assumed to be 21% of the total cell volume, and thus the thickness of the outer shell in this model is 2.7 µm. The surface or cortical tension coefficient for the interface between the cytoplasm and the surrounding fluid (plasma) is assumed to be 3.0×10 -5 N/m. The viscosity of the cytoplasm is higher than that of the plasma. The cell radius is assumed to be 6.5 µm, as in [START_REF] Khismatullin | A 3D numerical study of the effect of channel height on leukocyte deformation and adhesion in parallel-plate flow chambers[END_REF], which is larger than the commonly accepted mean radius of a neutrophil that is around 4.3 µm [START_REF] Needham | Rapid flow of passive neutrophils into a 4 µm pipet and measurement of cytoplasmic viscosity[END_REF]. In the second model shown in Figure 2b, called the compound-drop model, the outer layer is the same as above but the nucleus is modeled by a viscoelastic liquid. Thus, in addition to the outer interface between the shell and the surrounding fluid, there is an interface between the core fluid and the shell. The surface tension for both interfaces is assumed to be 3.0×10 -5 N/m and the viscosity of the viscoelastic liquid inside the nucleus is assumed to be higher than that of the shell. The values of the parameters for the two previous models are given in Table 1. The densities of the plasma and leukocyte are assumed to be equal, and thus the gravity plays no role in the dynamics.

Adhesion models

Under normal conditions, leukocytes are freely suspended in the blood stream. However, upon inflammation, the leukocytes that come close to the endothelial layer, as a result of the collisions with other blood cells and the expression or activation of adhesion molecules, can be captured. This latter process can be modeled by an adhesion force F b , which acts between the leukocytes and the endothelium. This force has a very short range and can result in the capture of leukocytes that come sufficiently close to the wall. We next describe two models that are used to represent the adhesion force.

Table 1. Parameters for the two-layered leukocyte models. Parameters marked with a star apply to the compounddrop model only. Appropriate references for the choice of these values have been included in brackets.

Adhesion potential model

In the first model, introduced by Sukumaran and Seifert [START_REF] Sukumaran | Influence of shear flow on vesicles near a wall[END_REF], the adhesion potential of the leukocyte with a plane wall at x = 0 is assumed to be

2 2 0 0 ( ) 2 d d W x w x x = -                   (1)
where w is the adhesion strength and 0 0.02

d R =
is a constant related to the cell radius R. Depending on its magnitude, the latter models both specific and nonspecific interactions.

Kinetic adhesion model

The second adhesion model is the kinetic model of Dembo et al. [START_REF] Dembo | The reaction-limited kinetics of membrane-to-surface adhesion and detachment[END_REF]. Here, the adhesion force between the cell membrane and the wall is assumed to arise due to the formation of bonds between the specific adhesion molecules on the cell membrane and the receptors on the substrate. A bond between the cell and the substrate is modeled as a Hookean spring, and the corresponding force f b due to one bond is ( )

b m mv f x l σ λ = -- (2
) where σ is the spring constant, x m is the distance between the cell surface and the wall, l mv is the unstressed microvillus length [START_REF] Khismatullin | Three-dimensional numerical simulation of a receptor-mediated leukocyte adhesion to surface: Effects of cell deformability and viscoelasticity[END_REF][START_REF] Bongrand | Cell Surface Dynamics: Concepts and Models[END_REF][START_REF] Hammer | Simulation of cell rolling and adhesion on surfaces in shear flow: General results and analysis of selectin-mediated neutrophil adhesion[END_REF][START_REF] Shao | Static and dynamic lengths of neutrophil microvilli[END_REF], and λ is the equilibrium or unstretched bond length. The total force per unit surface area F b between the cell and the substrate is then deduced from F N f =

(3) where N b is the bond density per unit area. The time evolution of the latter is given by a kinetic equation which balances the formation and dissociation of bonds:

0 0 ( )( ) b f l b r b r b N k N N N N k N t ∂ = - -- ∂ , (4)
where N l0 and N r0 are the initial ligands and receptors densities on the surface of the cell membrane.

The reverse and forward reaction rate coefficients k r and k f in (4) are given by: 100.0 Pa.s [START_REF] Drury | Hydrodynamics of micropipette aspiration[END_REF][START_REF] Tsai | Passive mechanical behavior of human neutrophils: Power-law fluid[END_REF] where k r0 and k f0 are the equilibrium values of the reverse and forward reaction rates, σ and σ ts are the spring and transition spring constants, k b is the Boltzmann constant, and T is the temperature. The values of these parameters are listed in Table 2.

2 0 ( )( ) exp 2 ts m mv r r b x l k k k T - -- =       σ σ λ (5) 2 0 ( ) exp 2 ts m mv f f b x l k k k T -- = -       σ λ (6)
The equilibrium bond density N b0 reaches the maximum value N b0m when the bond stretching x = x mλ -l mv is zero. For a given value of x m , the equilibrium bond density N b0 satisfies the following quadratic equation, which is deduced by setting the time dependent term in (4) equal to zero ( )( )

0 0 0 0 0 0 = - - - b r b r b l f N k N N N N k (7)
and which admits two solutions:

0 0 1,2 ( ) 2 r r l f k N N D k N + + ± = ( 8
)
where D is the discriminant of the quadratic equation, i.e.

2 0 0 0 0 ( ) 4 r r l r l f k D N N N N k = + + - (9)
The larger solution N 1 leading to a value larger than the number of available bonds is not physically meaningful. The only physical solution of ( 7) is thus N 2 , leading to the equilibrium bond density:

0 0 0 ( ) 2 r r l f b k N N D k N + + - = (10)
For the parameter values listed in Table 2, the maximum value of the equilibrium bond density, which is deduced by setting k f = k f0 and k r =k r0 , is N b0m =1.5×10 14 m -2 . Both the equilibrium bond density N b0 and the adhesion force vary with the distance between the cell membrane and the wall. In Figure 3a . Notice that when x' is approximately 48 the force reaches the maximum possible value. In terms of the dimensional variables, the force is maximum when |x m -λ -l mv | = 14.9 nm, and the adhesion force is approximately zero when |x m -λ -l mv | is greater than 27.9 nm. Therefore, the force is present only over the narrow distance range of ~30 nm. Consequently, as discussed later, a cell can be captured only if it is located within the above range.

Parameters

We next consider the time evolution of N b when x m is suddenly increased from (λ + l mv ) at t=0, to 1.01(λ + l mv ), which causes the bonds to break and the bond density to evolve to a new equilibrium value (see Figure 3b). The bond density N b in this case decays with time from N b0m to a smaller value. Our objective here is to determine the characteristic time over which the bond density evolves. This is important because in a flow situation the distance between the cell wall and the endothelial monolayer changes with time. Clearly, if this characteristic time over which the new equilibrium value is reached is much smaller than the flow time scale, it is appropriate to assume in the fluid dynamics problem that the bond density is given by the steady solution of equation ( 4). The bond density as a function of time when at t=0, x m is suddenly increased from (λ + l mv ) to 1.01(λ + l mv ). The new equilibrium value is approximately reached at t=7x10 -6 s.

By integrating equation ( 4) analytically, it is easy to show that

) / exp( 1 ) / exp( ) ( 1 2 τ β τ β t t N N t N - - - - = (11)
where N 1 and N 2 are the two roots defined in [START_REF] Bell | Cell adhesion. Competition between nonspecific repulsion and specific bonding[END_REF] and D is given by (

D τ = 0 1 2 < - - = N N α α β 9), 1 , f k
, and α is the initial bond density which, for the case considered, is

m b N 0 .
From equation [START_REF] Sukumaran | Influence of shear flow on vesicles near a wall[END_REF] it is clear that as

∞ → t
, we obtain N(t) = N 2 as the solution for the new equilibrium state. Furthermore, it is interesting to notice that for all t > 0, we have 0 )) / exp( 1 (

) / exp( ) ( ) ( ' 2 2 1 < - - - - = τ β τ τ β t t N N t N . ( 12
)
We then define the characteristic time t c by setting N 2 =N'(0) t c + α, which after using (12) gives,

1 2 1 1 ( ) 1 ( ) ( ) c f N N t N k N τ - = = -α -α (13)
From the above expression, we deduce that while t c depends on 0 0 , , ,

f r r l k k N N , it is independent of , ts σ σ .
For the parameter values listed in Table 2, the characteristic time t c given by equation ( 13) is O(10 -4 ) s. If the shear rate is between 10 and 100 s -1 , the characteristic flow time will be between 0.01~0.1 s. Clearly, only when these two time scales are comparable we need to consider the time evolution of the bond density for computing the adhesion force. If this is not the case, which happens for the parameter values selected in this paper, it is appropriate to assume in the fluid dynamics problem that the bond density is equal to the equilibrium value for that particular bond length. The adhesion force, in this case, depends only on the distance between the cell surface and the wall, which is also the case for the potential model. The force range for the kinetic model is however much shorter.

GOVERNING EQUATIONS

Let us denote the domain containing the plasma and the leukocyte by Ω , and the domain boundary by Γ . The plasma is modeled as a Newtonian liquid while the nucleus is modeled as a rigid sphere P(t) or as a viscoelastic liquid. The governing equations for this system are

0 = ⋅ ∇ u (14) ρ ∂u ∂t + u⋅∇u       =ρg-∇p+∇⋅ c λ r A       +∇⋅ 2η s D ( ) +γκδ(φ)n + ) (φ δ b F u = u L on Γ, u = U + ω x r on ) t ( P ∂
The nucleus in the compound-drop model of leukocytes is assumed to be a viscoelastic liquid modeled by the Oldroyd-B model. The viscoelastic stresses, in terms of the configuration tensor A, is given by

T r 1 ( ) t ∂ + ⋅∇ = ⋅∇ +∇ ⋅ - - ∂ λ A u A A u u A A I ( 15
)
Here u is the velocity, p is the pressure, η s is the solvent viscosity, ρ is the density, D is the symmetric part of the velocity gradient tensor, c is the polymer concentration parameter, λ r is the relaxation time of the viscoelastic fluid, n is the outer normal, γ is the surface tension, κ is the surface curvature, φ is the level set function which is defined to be the distance from the interface and δ is the delta function . A code based on the finite element method [START_REF] Pillapakkam | A Level-Set Method for Computing Solutions to Viscoelastic Two-Phase Flow[END_REF][START_REF] Singh | Finite element simulation of the start-up problem for a viscoelastic problem in an eccentric cylinder geometry using third-order upwind scheme[END_REF][START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Singh | Direct numerical simulation of viscoelastic particulate flows[END_REF][START_REF] Singh | Fluid dynamics of Floating particles[END_REF][START_REF] Kadaksham | Dynamics of pressure driven flows of electrorheological suspensions subjected to spatially non-uniform electric fields[END_REF] is used for solving the deformation of a leukocyte in a pressure driven flow. Here, the governing equations are solved simultaneously both inside and outside the leukocyte. The cell interface is tracked using the level set method [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF][START_REF] Pillapakkam | A Level-Set Method for Computing Solutions to Viscoelastic Two-Phase Flow[END_REF] and the adhesion force between the wall and the leukocyte is treated as a body force.

RESULTS

We next describe the direct simulation results for the motion of leukocytes in pressure driven flows in three dimensional channels with rectangular cross sections. The distance h between the left and right walls of the channel is 50 µm, the width of the channel is 37 µm, and the radius of the leukocytes is 6.5 µm [START_REF] Kan | Hydrodynamics of a compound drop with application to leukocyte modeling[END_REF][START_REF] N'dri | Computational modeling of cell adhesion and movement using a continuum-kinetics approach[END_REF][START_REF] Khismatullin | Three-dimensional numerical simulation of a receptor-mediated leukocyte adhesion to surface: Effects of cell deformability and viscoelasticity[END_REF]. All lengths reported in this paper have been nondimensionalized with respect to the leukocyte radius. For most cases reported, the dimensionless parameters are Re=0.5, Ca=0.167, De= 0.1 and Adh=50.0. For some cases Re, Ca and Adh are varied.

Leukocyte motion in the absence of adhesion forces.

We first consider the case where the adhesion force between the wall and the leukocytes is not present, and thus the behavior of a leukocyte is similar to that of a drop in a pressure driven flow. In a pressure driven flow, since the shear stress is maximum at the channel wall and decreases linearly with increasing distance from the wall, the deformation is expected to be maximum when the leukocyte is close to the wall and to decrease with increasing distance from the wall. The stretching occurs at an angle of 45 o with the flow direction, and the stretched material elements rotate due to the flow vorticity.

Compound-drop model. Figure 4 shows the steady state leukocyte shapes for the compound-drop model at various distances from the channel wall. In all cases, both the outer layer and the nucleus of the leukocyte are deformed. Also, as noted above, the deformation of the leukocyte increases, and the angle between the principal axis of the deformed shape and the channel wall decreases, as the cell gets closer to the channel wall. In addition, the angle between the principal axis of the deformed nucleus and the wall is greater than that between the deformed cell and the wall. Moreover, the viscoelastic stresses inside the nucleus are maximal in the region which is closest to the wall. Drop-rigid-particle model. In the drop-rigid-particle model, the nucleus (a rigid sphere) remains spherical while the outer layer of the leukocyte deforms. The steady deformed shapes for the drop-rigid-particle model at various distances from the wall are shown in Figure 5. Figures 4 and5 show that the overall deformation is smaller for the drop-rigid-particle model. This implies that the drop-rigid-particle model is not appropriate for situations in which the cell, as well as the nucleus, deform significantly, such as during spreading or transmigration.

Figure 6 displays the steady shape of a leukocyte released at a distance of x=2 from the wall. It is stretched and its principal direction makes an angle of ~40 degrees with the channel wall. Our simulations show that the leukocyte is moving laterally, away from the channel wall, while moving in the flow direction. This is also the case for the compound-drop model described in the previous subsection. Lateral migration of leukocytes. Past studies of the motion of drops in the pressure driven flows have shown that a drop not only moves in the flow direction, but also drifts laterally [START_REF] Segré | Behavior of macroscopic rigid spheres in Poiseuille flow. Part 2. Experimental results and interpretation[END_REF][START_REF] Goldsmith | The flow of model particles and blood cells and its relation to thrombogenesis[END_REF]. The velocity of the lateral drift, however, is very small compared to the velocity in the flow direction. Therefore, regardless of the initial position of the drop, it eventually, after a relatively large displacement in the flow direction, reaches a lateral equilibrium position, which depends on the ratio of the viscosity of the fluid inside the drop and that of the surrounding fluid.

For the parameter values considered in this study, the Reynolds number is 0.5, and the ratio of the cytoplasm and plasma viscosities is 100. The lateral migration of leukocytes, released at different distances from the wall of the channel, is studied in a periodic domain. The initial positions ranged between x= 1.475 and x=3.375. The channel wall is at x=0 and the centerline of the channel is at x=3.85. Some of the numerically obtained trajectories of leukocytes' centers of mass are displayed in figure 7.

Our simulations show that the lateral equilibrium position of a leukocyte is around x=2.5, which is approximately at a distance of 0.32 h from the channel wall, where h is the channel height. This implies that a leukocyte cannot be captured at a wall unless it is brought closer to the wall by some other mechanism, e.g., by collisions with other blood cells.

Leukocyte motion in the presence of adhesion forces.

We next describe the motion of a leukocyte released near the channel wall at a distance such that the adhesion force is present. Our objective here is to study its motion and deformation under the combined action of the adhesive and hydrodynamic forces. We consider both models of the adhesive force, i.e., the potential model ( 1) and the kinetic model described in section 2.4.

In the presence of a pressure driven flow, the leukocyte experiences a hydrodynamic lift force whose magnitude depends on many factors including the flow velocity. If the lift force is larger than the adhesion force, the leukocyte moves away from the wall. Obviously, when the flow velocity is smaller than a critical value such that the lift force is smaller than the adhesive force, the leukocyte is expected to remain attached to the wall. Adhesion potential model. In a pressure driven flow, in addition to the adhesion force the leukocyte also experiences a hydrodynamic force. The component of the hydrodynamic force in the flow direction, or drag, moves it in the flow direction and the component normal to the wall, or lift, moves it away from the wall, thus countering the adhesion force acting towards the wall. The adhesion potential strength and the initial leukocyte position of x=2.0 were selected so that the leukocyte is attracted to the wall. We first consider the case of the drop-rigid-particle model.

The flow causes the leukocyte to stretch along the shear direction, while the adhesive force pulls it towards the wall (see Figure 8). The adhesive force in the potential model is long-ranged, which can be seen by the fact that a leukocyte located at a short distance from the wall is pulled towards the wall. This feature of the potential model is perhaps unrealistic since the adhesive forces are not long-ranged. We remind the reader that the force for the kinetic model is very short-ranged. As the trailing end of the deformed leukocyte is closer to the wall, it experiences a larger attractive potential force, and thus is pulled further towards the wall. Consequently, the trailing end of the leukocyte begins to flatten, and the contact area between the leukocyte and the wall starts increasing. The contact area continues to increase until the surface tension force becomes equal to the deforming forces.

For the compound-drop model, the nucleus of the leukocyte also deforms (Figure 9). The other features of the overall deformation are similar. Here, the initial position of the center of the leukocyte is x=0.5. The flow causes the leukocyte to move upwards, and the adhesive and shear forces cause the surface of the leukocyte near the wall to flatten. Furthermore, it is interesting to note that the gap between the leukocyte surface and the wall near the leading edge is larger than that near the trailing edge. Notice that this flattening of the leukocyte cannot be seen in the top view where it appears to be only elongated [START_REF] Khismatullin | Three-dimensional numerical simulation of a receptor-mediated leukocyte adhesion to surface: Effects of cell deformability and viscoelasticity[END_REF][START_REF] Khismatullin | A 3D numerical study of the effect of channel height on leukocyte deformation and adhesion in parallel-plate flow chambers[END_REF]. This is an important point because in experiments we generally have access only to the top view. Kinetic adhesion model. As noted above, the adhesion force for the kinetic model is non-zero for a very small range of distances between the leukocyte surface and the wall. Therefore, for the results presented in this section the initial position of the leukocyte is selected so that the force is present. Figure 10 shows that at the beginning the leukocyte rolls near the wall surface. The Reynolds number for this case is 0.5. As the leukocyte deforms, the lift force it experiences increases and in this case the leukocyte moves away from the wall as the adhesion force is not sufficiently large to keep it attached. After the leukocyte moves out of the range of the kinetic adhesion force, which as noted above corresponds to the relatively short distance of about 30 nm, the adhesion force becomes zero. The leukocyte then migrates to a lateral distance at which the lateral hydrodynamic force is zero, as discussed above. In order to ensure that the leukocyte is captured at the wall, we next consider the case where the adhesion force is 100 times larger and the adhesion parameter is 5000. In this case, the cell remained attached to the wall (Figure 11). After the steady state is reached and the deformation no longer increases, the lower area of the cell is flattened, and the leukocyte rolls and slides on the adhesive wall. In addition, the front edge of the leukocyte is farther away from the wall and the trailing end (located near the wall) is flattened. The flattening increases the surface area over which the adhesion force acts, which is especially critical in the kinetic model, as the distance range over which the force acts is rather small [START_REF] Khismatullin | A 3D numerical study of the effect of channel height on leukocyte deformation and adhesion in parallel-plate flow chambers[END_REF][START_REF] Cantat | Dynamics and similarity laws of adhering vesicles in haptotaxis[END_REF].

DISCUSSION AND CONCLUSIONS

We have performed three-dimensional direct numerical simulations of the motion and deformation of a leukocyte in channels. The surface of the leukocyte is tracked using the level set method and in the case of the solid nucleus model the rigid motion enforcement inside the nucleus is performed by using the distributed Lagrange multiplier method. The leukocyte is represented using two different models: the compound-drop model and the drop-rigid-particle model. Preliminary results using this approach were presented in [START_REF] Jin | Direct simulation of the migration of leucocytes in pressure driven flow[END_REF].

Our simulations show that the leukocyte deformation is larger when it is closer to the wall, where the shear rate is larger. The deformation is smaller for the drop-rigid particle model of the leukocyte. For the parameter values selected, which are typical for such flows, the leukocyte moved away from the wall to a stable location positioned between the center of the channel and the wall. This implies that the adhesion force is necessary for the leukocyte to remain attached to the wall.

In the presence of a pressure driven flow, the leukocyte experienced not only a lift force whose magnitude depends on many parameters including the flow velocity, but also a hydrodynamic torque. Both cause the trailing end of the leukocyte to move away from the wall. Clearly, when the bulk flow velocity is sufficiently small the combined effect of these hydrodynamic forces is smaller than the adhesive force, and the leukocyte remained attached to the wall. However, when the flow velocity was larger and closer to typical values, the hydrodynamic forces overcame the adhesive force and the leukocyte moved away from the wall. Our simulation results showed that for the kinetic model parameter values selected in section 2.4 and typical flow rates, the leukocyte was found to move away from the wall. This is somewhat problematic since in experiments leukocytes do get attached to the endothelium monolayer.

For the two adhesion models used in this study the trailing end of the leukocyte flattened. The flattening increases the contact surface area over which the adhesion force acts, and thus the adhesion force acting on the leukocyte. The flattening, therefore, plays an important role in increasing the adhesion force so that the latter does not get canceled by the hydrodynamic lift force and torque. The flattening of the leukocyte is particularly crucial for the kinetic model since the distance range over which the corresponding adhesive force acts is approximately 30 nm.

Figure 2 .

 2 Figure 2. Two layered models of leukocytes. (a) drop-rigid-particle model, (b) compound-drop model.

  the adhesion force and bond density are shown as a function of the dimensionless

Figure 3 .

 3 Figure 3. Kinetic adhesion model. (a)The adhesion force F b and the equilibrium bond density N b0 are shown as functions of the dimensionless bond stretching x ' . (b)The bond density as a function of time when at t=0, x m is suddenly increased from (λ + l mv ) to 1.01(λ + l mv ). The new equilibrium value is approximately reached at t=7x10 -6 s.

  . The polymer contribution to the viscosity is η p = cη S , and the zero shear viscosity of the fluid is η o = η s + η p .The equations for the fluid and leukocyte motions are nondimensionalized by taking the characteristic length, velocity, time and stress scales as R, U, R/U and η 0 U/R, respectively. Here R is the leukocyte radius, 0 η is the plasma viscosity, and U is the centerline velocity of the undisturbed pressure driven flow. The nondimensional parameters are: Reynolds number

Figure 4 .

 4 Leukocyte deformation in the case of the compound-drop model at different distances from the wall. The shape on the xz-mid plane through the leukocyte center are shown. (a) x=1.9, (b) x=2.425, (c) x=2.9, and (d) x=3.375 at time t=4.6. Notice that both the leukocyte and the nucleus are deformed. The Reynolds number is 0.5.

Figure 5 .

 5 Leukocyte deformation in the case of the drop-rigid-particle model versus the distance from the wall. (a) x=1.9, (b) x=2.425, (c) x=2.9, and (d) x=3.375 at time t=4.6. The deformation increases with decreasing gap between the leukocyte and the wall. Since the nucleus is rigid, it does not deform. The Reynolds number is 0.5.

Figure 6 .

 6 Figure 6. Leukocyte deformation in the case of the drop-rigid-particle model at (a) t=0.82 and (b) t=5.25. The Reynolds number is 0.5.

Figure 7 .

 7 Figure 7. Trajectories of leukocytes released from different lateral positions within the channel. The Reynolds number is 0.5.

Figure 8 .

 8 Deformation of a leukocyte in the case of the drop-rigid-particle model due to shear forces while it moves towards the wall under the action of the adhesive force at (a) t=0.82 and (b) t=2.46. The Reynolds number is 0.5, Ca=0.167, De= 0.1 and Adh=50.0.

Figure 9 .

 9 Figure 9. Leukocyte deformation for the compound-drop model released close to the wall at (a) t=0.82, and (b) t=2.46, and (c) t=3.94. The top view is also shown. The Reynolds number is 0.5, Ca=0.167, De= 0.1 and Adh=50.0.

Figure 10 .Figure 11 .

 1011 Figure 10. Leukocyte deformation in the case of the compound-drop model and the adhesion force given by the kinetic model at (a) t=0.6, (b) t=3, and (c) t=6, for a Reynolds number of Re=0.5, Ca=0.167, De= 0.1 and Adh=50.0.

Table 2 .

 2 Parameter values for the kinetic adhesion model. Appropriate references have been included in brackets.

	Transition spring constant,	ts σ , N/m	-0.005 -0.0095[8]	10 -4
	Unstressed microvillus length, l mv , m	3×10 -7 [19, 21]	3×10 -7
	Equilibrium bond length, λ , m	10 -8 ~10 -7 [19]	10 -8
	Thermal energy, k b T, N•m		3.8 ~ 4.3×10 -21	4.28×10 -21


Martin J Gander

Laurence Halpern

Johann Rannou

Juliet Ryan

A DIRECT TIME PARALLEL SOLVER BY DIAGONALIZATION FOR THE WAVE EQUATION

With the advent of very large scale parallel computers, it has become more and more important to also use the time direction for parallelization when solving evolution problems. While there are many successful algorithms for diffusive problems, only some of them are also effective for hyperbolic problems. We present here a mathematical analysis of a new method based on the diagonalization of the time stepping matrix proposed by Maday and Rønquist in 2007. Like many time-parallelization methods, at first this does not seem to be a very promising approach: the matrix is essentially triangular, or, for equidistant time steps, actually a Jordan block, and thus not diagonalizable. If one chooses however different time steps, diagonalization is possible, and one has to trade off between the accuracy due to necessarily having different time steps, and numerical errors in the diagonalization process of these almost non-diagonalizable matrices. We present for the first time such a diagonalization technique for the Newmark scheme for solving wave equations, and derive a mathematically rigorous optimization strategy for the choice of the parameters in the special case when the Newmark scheme becomes Crank-Nicolson. Our analysis shows that small to medium scale time parallelization is possible with this approach. We illustrate our results with numerical experiments for model wave equations in various dimensions, and also an industrial test case for the elasticity equations with variable coefficients.

1. Introduction. Using the time direction in evolution problems for parallelization is an active field of research; see [START_REF]50 years of time parallel time integration[END_REF] for an overview. Most of the methods developed for this purpose are iterative, see for example the parareal algorithm [START_REF] Lions | A parareal in time discretization of PDE's[END_REF], whose convergence was analyzed in [START_REF] Gander | Analysis of the parareal time-parallel timeintegration method[END_REF] for linear problems. It was also shown in [START_REF] Gander | Analysis of the parareal time-parallel timeintegration method[END_REF] that for the transport equation no speedup is possible with the parareal algorithm; for further analyses in the hyperbolic case, see [START_REF] Gander | Analysis of the parareal algorithm applied to hyperbolic problems using characteristics[END_REF]. A sharp convergence estimate of the parareal algorithm for non-linear problems can be found in [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF]. Later, a variation of the parareal algorithm using spectral deferred correction [START_REF] Michael | A hybrid parareal spectral deferred corrections method[END_REF] led then to the PFASST algorithm [START_REF] Emmett | Toward an efficient parallel in time method for partial differential equations[END_REF], which is a multilevel method. Space-time multigrid methods were also developed (see [START_REF] Horton | A space-time multigrid method for parabolic partial differential equations[END_REF] and references therein), and a new such method using only standard components can be found in [START_REF] Gander | Analysis of a new space-time parallel multigrid algorithm for parabolic problems[END_REF] with excellent strong and weak scaling properties for parabolic problems. There is also the non-invasive MGRIT algorithm [START_REF] Friedhoff | A multigrid-in-time algorithm for solving evolution equations in parallel[END_REF][START_REF] Falgout | Parallel time integration with multigrid[END_REF], and it was shown in [START_REF] Gander | Multigrid interpretations of the parareal algorithm leading to an overlapping variant and MGRIT[END_REF] that MGRIT is equivalent to an overlapping parareal algorithm, which led to a detailed non-linear convergence analysis for MGRIT [START_REF] Gander | Multigrid interpretations of the parareal algorithm leading to an overlapping variant and MGRIT[END_REF]. Early versions of these space-time multigrid methods were based on waveform relaxation [START_REF] Lubich | Multi-grid dynamic iteration for parabolic equations[END_REF][START_REF] Vandewalle | Space-time concurrent multigrid waveform relaxation[END_REF], and there are also very successful Schwarz waveform relaxation methods for the time parallel solution of evolution partial differential equations [START_REF] Gander | Space time continuous analysis of waveform relaxation for the heat equation[END_REF][START_REF] Giladi | Space time domain decomposition for parabolic problems[END_REF][START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF][START_REF] Gander | Overlapping Schwarz waveform relaxation for convectiondominated nonlinear conservation laws[END_REF][START_REF]Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Bennequin | A homographic best approximation problem with application to optimized Schwarz waveform relaxation[END_REF], and the more recent Dirichlet-Neumann and Neumann-Neumann waveform relaxation methods [START_REF] Mandal | A time-dependent Dirichlet-Neumann method for the heat equation[END_REF][START_REF] Kwok | Neumann-Neumann waveform relaxation for the time-dependent heat equation[END_REF][START_REF] Gander | Dirichlet-Neumann and Neumann-Neumann waveform relaxation algorithms for parabolic problems[END_REF]. These are among the very few space-time iterative methods that are effective for hyperbolic problems [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF][START_REF] Gander | Absorbing boundary conditions for the wave equation and parallel computing[END_REF][START_REF] Qaddouri | Optimized Schwarz methods with an overset grid for the shallow-water equations: preliminary results[END_REF][START_REF] Dirichlet-Neumann | Neumann waveform relaxation for the wave equation[END_REF], since they use the finite speed of propagation for good space-time decompositions; see also the related tent-pitching approach [START_REF] Gopalakrishnan | Mapped tent pitching schemes for hyperbolic systems[END_REF] which removes the need for iteration completely by following in the decomposition the characteristics. A different approach for hyperbolic problems are the Krylov parareal methods [START_REF] Farhat | Timeparallel implicit integrators for the near-real-time prediction of linear structural dynamic responses[END_REF][START_REF] Gander | Analysis of a modified parareal algorithm for second-order ordinary differential equations[END_REF][START_REF]Analysis of a Krylov subspace enhanced parareal algorithm for linear problems[END_REF][START_REF] Ruprecht | Explicit parallel-in-time integration of a linear acousticadvection system[END_REF], which however have an important overhead due to orthogonalization.

As an alternative, one can use direct space-time parallel solvers, like RIDC [START_REF] Christlieb | Parallel high-order integrators[END_REF], which is ideal for multicore architectures with 4 to 8 cores, and generates high order solutions in time in parallel at the runtime cost of an Euler time stepper. For linear hyperbolic problems, there is also the ParaExp algorithm [START_REF] Gander | Paraexp: A parallel integrator for linear initial-value problems[END_REF] which is based on a fully overlapping decomposition and Krylov techniques. A further direct approach based on the diagonalization of the time stepping matrix was introduced in [START_REF]Parallelization in time through tensor-product space-time solvers[END_REF]. This approach contains however parameters which must be carefully selected for the method to be successful. A first estimate for these parameters for a Backward Euler discretization of the heat equation was obtained in the short manuscript [START_REF] Gander | A direct solver for time parallelization[END_REF]. The idea of diagonalization can already be found in [START_REF] Bickart | An efficient solution process for implicit Runge-Kutta methods[END_REF][START_REF] Butcher | On the implementation of implicit Runge-Kutta methods[END_REF] at the level of implicit Runge-Kutta methods, but not with parallelization in mind. Wave equations are however more delicate, and require more refined techniques leading to more accurate estimates for the parameters.

We present a new such direct time parallel method based on diagonalization for wave equations discretized by the Newmark scheme. Our analysis for the special case of Crank-Nicolson shows that with this approach not an arbitrarily large number of processors can be used to parallelize in time, only about 3-4 times more than for RIDC for example, but as for RIDC, this number becomes a multiplicative factor for the processors used for parallelization in space: for example if 1000 processors were used for parallelization in space, and diagonalization permits the use of 20 processors in time, overall 20000 processors could then be used to solve the wave equation in parallel, instead of only 1000 without parallelization by diagonalization. We consider as our model partial differential equation (PDE) the second order wave equation, (

∂ tt u -∆u = 0 in Ω × (0, T ), u = 0 on ∂Ω, (u, ∂ t u)(•, 0) = (f, g) in Ω. 1.1)
A popular implicit time integrator for the wave equation (1.1) is the Newmark scheme [START_REF] Newmark | A method of computation for structural dynamics[END_REF]: let ∆ h be a space discretization of the Laplacian ∆, and we thus want to discretize üh -∆ h u h = 0 in time, where üh := ∂ tt u h . In the Newmark scheme, the additional unknown uh := ∂ t u h is introduced, and one then approximates both u h and uh on the time partition 0 = t

0 < t 1 < t 2 < • • • < t N = T , k n := t n -t n-1
, by a time stepping rule involving the two parameters β and γ, (1.2)

u n+1 h = u n h + k n+1 un h + k 2 n+1 (( 1 2 -β)ü n h + β ün+1 h ), un+1 h = un h + k n+1 ((1 -γ)ü n h + γ ün+1 h ), ün h -∆ h u n h = 0,
with given initial data u 0 h and u0 h . Let I t be the N ×N identity matrix associated with the time domain and I x be the J × J identity matrix associated with the spatial domain. Setting u := (u 1 h , . . . , u N h ), and similarly for u and ü, and defining the matrices

B 1 =      1 k1 -1 k2 1 k2 0 0 . . . . . . -1 k N 1 k N      , B 2 =      0 1 0 0 0 . . . . . . 1 0      , C 1 =      βk 1 ( 1 2 -β)k 2 βk 2 0 0 . . . . . . ( 1 2 -β)k N βk N      , C =      γ (1 -γ) γ 0 0 . . . . . . (1 -γ) γ     
, the scheme (1.2) can be written in the compact form

(1.3) (B 1 ⊗ I x )u = (B 2 ⊗ I x ) u + (C 1 ⊗ I x )ü + F, (B 1 ⊗ I x ) u = (C ⊗ I x )ü + G, ü -(I t ⊗ ∆ h )u = 0,
with the right hand sides

F = ( 1 k 1 u 0 h + u0 h + ( 1 2 -β)k 1 ∆ h u 0 h , 0, • • • , 0), G = ( 1 k 1 u0 h + (1 -γ)∆ h u 0 h , 0, • • • , 0).
Solving the second equation in (1.3) for u using that for the Kronecker product (B

1 ⊗ I x )(B 2 ⊗ I x ) = (B 1 B 2 ⊗ I x ), we get u = (B -1 1 ⊗ I x )((C ⊗ I x )ü + G)
, which we insert into the first equation to obtain

(B 1 ⊗ I x )u = ((B 2 B -1 1 C + C 1 ) ⊗ I x )ü + (B 2 B -1 1 ⊗ I x )G + F.
The matrix B 3 := B 2 B -1 1 C + C 1 is a lower triangular matrix whose diagonal is that of C 1 . Hence B 3 is invertible, and we can solve for ü to obtain ü = (B -1

3 ⊗ I x ) (B 1 ⊗ I x )u -(B 2 B -1 1 ⊗ I x )G -F .
Inserting this into the last equation in (1.3), we find

(1.4) (B ⊗ I x -I t ⊗ ∆ h )u = f , with (1.5) B = B -1 3 B 1 , f = (B -1 3 ⊗ I x )((B 2 B -1 1 ⊗ I x )G + F ).
Equation (1.4), is the matrix tensor-product form of the Newmark scheme for the wave equation in mixed form. It is equivalent to (1.2), and since the matrices are well-conditioned, it gives numerically the same solution. The velocity vector can be recovered using the equation

u = (B -1 1 ⊗ I x )((C ⊗ ∆ h )u + G).
The key idea to solve a space-time problem written in tensor-product form (1.4) in a time parallel fashion from [START_REF]Parallelization in time through tensor-product space-time solvers[END_REF] is to diagonalize B. Since all matrices are lower triangular, B is lower triangular, and the diagonal of B equals (

1 βk 2 1 , • • • , 1 βk 2 N ). If all the time steps are different, then B is indeed diagonalizable, B = SDS -1
, where the diagonal matrix D contains the eigenvalues 1 βk 2 n on the diagonal, and (1.4) can be written as

(1.6) (S ⊗ I x )(D ⊗ I x -I t ⊗ ∆ h )(S -1 ⊗ I x )u = f .
One can then solve (1.6) parallel in time performing the following 3 steps:

(1.7) (a) g = (S -1 ⊗ I x )f , (b) ( 1 βk 2 n -∆ h )w n = g n , 1 ≤ n ≤ N, (c) u = (S ⊗ I x )w.
The main work, namely the N linear systems in space in step (b), can all be solved in parallel using N processors, or if the spatial problems were already parallelized using M processors, one could then overall use N ×M processors. We call this approach time parallelization by diagonalization. It was first introduced and analyzed numerically for wave propagation in [START_REF] Rannou | PRF transition statique-dynamique dans les structures composites[END_REF][START_REF]PRF transition statique-dynamique dans les structures composites, rapport technique[END_REF], motivated by the original invention in [START_REF] Maday | Rønquist, Fast tensor product solvers. Part II: Spectral discretization in space and time[END_REF][START_REF]Parallelization in time through tensor-product space-time solvers[END_REF] for diffusive problems the authors called tensor-product space-time solvers. We will analyze here specifically the Crank-Nicolson scheme, which corresponds to β = 1 4 and γ = 1 2 in the Newmark scheme (1.2), and leads to the unconditionally stable time-integration operator of maximum accuracy most often used in practice. To see in what sense the parameter choice β = 1 4 and γ = 1 2 corresponds to Crank-Nicolson, we first replace the last equation in (1.2) into the first two for these values of the parameters,

u n+1 h = u n h + k n+1 un h + k 2 n+1 4 (∆ h u n h + ∆ h u n+1 h ), un+1 h = un h + kn+1 2 (∆ h u n h + ∆ h u n+1 h ).
Now from the second equation, we have

1 4 (∆ h u n h + ∆ h u n+1 h ) = ( un+1 h -un h )/(2k n+1
), which when inserted into the right hand side of the first equation gives

(1.8) u n+1 h = u n h + kn+1 2 ( un h + un+1 h ), un+1 h = un h + kn+1 2 (∆ h u n h + ∆ h u n+1 h ),
which is clearly Crank-Nicolson applied to the mixed formulation of the wave equation, which with v h := uh is given by

uh = v h vh = ∆ h u h .
Crank-Nicolson is second order accurate in time, and unconditionally stable. With a similar computation as for the Newmark scheme, we also obtain the tensor-product form (1.4) for the Crank-Nicolson scheme with (1.9)

B = (C -1 B 1 ) 2 , f = (C -1 ⊗ I x )((B 1 C -1 ⊗ I x )F + G), β = 1 4 , γ = 1 2
, and the same G, but with F = (

1 k1 u 0 h + 1 2 u0 h , 0, • • • , 0). 2.
Analysis of time parallelization by diagonalization. Using the time parallel solver (1.7) based on diagonalization requires some care: first, the time stepping matrix B is only diagonalizable if the time steps are all different, and this can lead to a larger discretization error compared to using equidistant time steps. Second, the condition number of the eigenvector matrix S increases exponentially with the number of time steps N to be done in parallel, which leads to inaccurate results in steps (a) and (c) of (1.7) because of roundoff error. One therefore needs to carefully and accurately estimate these two errors, to determine how many time steps can indeed be performed in parallel using (1.7) without loosing accuracy.

In order to obtain such estimates in the wave equation case, we use a Fourier transform in space with Fourier variable ξ, and obtain ∂ 2 û ∂t 2 +|ξ| 2 û = 0, which indicates that we need to study for fixed T and a > 0 the algorithm (1.7) applied to the ordinary differential equation (ODE) (2.1) ü + a 2 u = 0, t ∈ (0, T ), with initial conditions u(0) = f and u(0) = g. The Crank-Nicolson scheme for (2.1) is Geometric time stepping mesh for q = 1.1, a typical value of the parameter we will determine with our analysis, compared to an equidistant time stepping mesh with initial conditions u 0 = f and u0 = g. Let U := (u,1 a u) be the solution of (2.1) including the derivative term which will be important for our estimates, and U n := (u n , 1 a un ) be the discrete approximation defined by (2.2). We can then rewrite the continuous and discretized problems as

(2.2) u n = u n-1 + 1 2 k n ( un + un-1 ), un = un-1 -a 2 2 k n (u n + u n-1
(2.3) d t U + aJU = 0, (I+ ak n 2 J)U n = (I- ak n 2 J)U n-1 , J := 0 -1 1 0 , with initial condition U (0) = U 0 = (f, g a ).
For a given N , we now consider time steps given by a geometric partition 1 T q := (k 1 , • • • , k N ), k n (q) := q n-1 k 1 as it was suggested in [START_REF] Maday | Rønquist, Fast tensor product solvers. Part II: Spectral discretization in space and time[END_REF]. The constraint N n=1 k n = T determines uniquely the value of the geometric time steps, (2.4) k n = q n N j=1 q j T.

We call U n (T q ) the corresponding sequence of solutions. The equidistant timestep k = T N is obtained for q = 1. An example of such a geometric time mesh is given in Figure 2.1 for a typical value q = 1.1 our analysis will determine, and for comparison also an equidistant time mesh is shown.

2.1. Estimation of the error due to the geometric time stepping. We now present an accurate estimate for the change in the truncation error caused by using a geometric time stepping with (2.5) q = 1 + ε, ε a tunable parameter, compared to using equal time stepping, by studying the distance between U N (T q ) and U N (T 1 ). Theorem 2.1 (Asymptotic truncation error estimate). For fixed a, T and N , we obtain for ε small the error estimate

(2.6) U N (T 1+ε ) -U N (T 1 ) = φ( aT 2N , N )ε 2 U 0 + O(ε 3 ),
where φ(y, N ) := N (N 2 -1) 6 y 3 (1+y 2 ) 2 . Proof. We first diagonalize the two systems (2.3) in the complex plane: let

P := 1 1 i -i , D := -1 0 0 1 ,
which implies that P -1 JP = iD and P -1 (I -αJ)P = I -iαD. Defining

µ(t) := e -iat , λ n := 1 -i akn 2 1 + i akn 2
, and

µ n := n p=1 λ p ,
we obtain the diagonalized continuous and discrete solutions

U (t) = P μ(t) 0 0 µ(t) P -1 U 0 , U n = P μn 0 0 µ n P -1 U 0 .
For a given N , we now estimate the change induced by the non equidistant time steps.

For any partition

T = (k 1 , • • • , k N ) of the interval (0, T ), we compare the solution U N (T ) to the solution U N ( T ) obtained by the equidistant mesh T = ( k, • • • , k) with k = T /N . Defining µ(T ) := N n=1 1 -i akn 2 1 + i akn 2 ,
we find that

(2.7) U N (T ) -U N ( T ) = P µ(T ) -µ( T ) 0 0 µ(T ) -µ( T ) P -1 U 0 .
Taking norms, and noting that P V 2 = 2 V 2 and P -1

V 2 = 1 2 V 2 , we get (2.8) U N (T ) -U N ( T ) = |µ(T ) -µ( T )| U 0 , U N ( T ) -U (T ) = |µ( T ) -e -iaT | U 0 .
The second formula is the truncation error at the end of the time interval for equidistant timesteps. We next need to estimate the error in µ, for which we use the following result. Lemma 2.2. For fixed a, T and N , we have

(2.9) |µ(T ) -µ( T )| = a 3 k 4(1 + a 2 k2 4 ) 2 T -T 2 + O(||T -T || 3 ).
Proof. We write µ(T ) as

µ(T ) = N n=1 λ(k n ), λ(k) = 1 -i ak 2 1 + i ak 2 .
By the second order Taylor-Young formula, we have

µ(T ) = µ( T ) + µ ( T ) • (T -T ) + 1 2 µ ( T ) • (T -T ) • (T -T ) + O(||T -T || 3 ).
We now compute the derivative,

(2.10) ∂µ ∂k n ( T ) = µ( T ) λ (k n ) λ(k n ) = µ( T ) -ia 1 + a 2 k2 4 ,
and therefore obtain for the linear term

µ ( T ) • (T -T ) = µ( T ) -ia 1 + a 2 k2 4 n (k n -k) = 0.
For the quadratic term, we differentiate (2.10) again with respect to k n , (2.11)

∂ 2 µ ∂k 2 n ( T ) = µ( T )a 2 -1 + i a k 2 (1 + a 2 k2 4 ) 2 , ∂ 2 µ ∂k n ∂k j ( T ) = µ( T )a 2 -1 (1 + a 2 k2 4 ) 2 ,
and therefore obtain

µ ( T ) • (T -T ) • (T -T ) = -a 2 µ( T ) n,j (k n -k)(k j -k) -i a k 2 n (k n -k) 2 (1 + a 2 k2 4 ) 2
, which gives

µ(T ) = µ( T ) - a 2 µ( T ) 2 n,j (k n -k)(k j -k) -i a k 2 n (k n -k) 2 (1 + a 2 k2 4 ) 2 + O(||T -T || 3 ).
The first sum in the numerator is zero, and hence

µ(T ) -µ( T ) = iµ( T )a 3 k n (k n -k) 2 4(1 + a 2 k2 4 ) 2 + O(||T -T || 3 ),
which gives (2.9) when taking the modulus, since the quadratic term is purely imaginary and |µ( T )| = 1. It remains to estimate T -T 2 , which is done in Lemma 2.3. The difference between the geometric and equal time partition satisfies the estimate

(2.12) T -T 2 = N n=1 (k n -k) 2 = k2 N (N 2 -1) 12 ε 2 + O(ε 3 ).
Proof. The time step k n in (2.4) has for ε small the expansion

k n = k(1 + α n ε + β n ε 2 + o(ε 2 )), with α n = n -N +1 2 and β n = n(n -N -2) + (N +1)(N +5) 6 . These coefficients satisfy the relations n α n = n β n = 0, n α 2 n = N (N 2 -1)

12

, which leads to (2.12).

We can now finish the proof of Theorem 2.1 by inserting (2.4) into (2.8) to obtain

U N (T ) -U N ( T ) = a 3 k 4(1 + a 2 k2 4 ) 2 k2 N (N 2 -1) 12 ε 2 U 0 + O(ε 3 ).
2.2. Estimation of the roundoff error due to diagonalization. The spacetime system (1.4) becomes in the scalar case (B + a 2 I)u = f , and to diagonalize this system, we need to diagonalize the matrix B, which will involve unipotent lower triangular Toeplitz matrices.

Definition 2.4. A unipotent lower triangular Toeplitz matrix of size N is of the form

(2.13) T (x 1 , . . . , x N -1 ) =       1 x 1 . . . . . . . . . . . . x N -1 . . . x 1 1      
.

The next Lemma gives a closed form eigendecomposition of the time stepping matrix B, and is proved using the techniques of q-hypergeometric series (see the book of Gasper and Rahman, Basic Hypergeometric series [START_REF] Gasper | Basic hypergeometric series[END_REF]). Lemma 2.5 (Eigendecomposition of B).

If k n = q n-1 k 1 with q = 1, then B has the eigendecomposition B = V DV -1 , with D = diag( 4 k 2 n
), and V and its inverse are unipotent lower triangular Toeplitz matrices given by

V = T (p 1 , . . . , p N -1 ), with p n := n j=1 1 + q j 1 -q j , (2.14)
V -1 = T (q 1 , . . . , q N -1 ), with q n := q -n n j=1 An eigenvector of C -1 B 1 is such that (B 1 -2 kn C)X (n) = 0, and the matrix B 1 -2 kn C is lower bidiagonal, which leads to a recursive formula for the coefficients X (n) of the eigenvectors,

1 + q -j+2 1 -q -j . ( 2
X (n) 1 = X (n) 2 = • • • X (n) n-1 = 0, X (n) n+1 = kn+kn+1 kn-kn+1 X (n) n , • • • , X (n) N = kn+k N kn-k N X (n) N -1 .
Since the normalization is arbitrary, we can choose

X (n) n = 1, which gives X (n) j =      0 for j < n, j l=n+1 kn+k l kn-k l for j > n.
The matrix of eigenvectors V := (X (1) ,

• • • , X (N )
) is lower triangular and unipotent, i.e. the diagonal elements equal 1. In case of the geometric mesh, it has the particular structure

V ij = X (j) i = i l=j+1 k j + k l k j -k l = i l=j+1 q j + q l q j -q l = i l=j+1 1 + q k-j 1 -q k-j = i-j l=1 1 + q l 1 -q l ,
which shows that V ij is a function of i-j only, and therefore V is a unipotent Toeplitz matrix. Consider now the inverse of V . First, it is easy to see that it is also unipotent Toeplitz. To establish (2.15) is equivalent to prove that

(2.16) for 1 ≤ n ≤ N -1, n j=0
p n-j q j = 0, with the convention that p 0 = q 0 = 1.

To show this result, we need to define Heine's q-hypergeometric series (Heine 1847),

(2.17)

2 ϕ 1 (a 1 , a 2 ; b; q; x) := ∞ n=0 (a 1 ; q) n (a 2 ; q) n (q; q) n (b; q) n x n ,
where we assume that b is not an integer negative power of q, and the q-rising factorial is

(a; q) k := k-1 i=0 (1 -q i a), (a; q) ∞ := ∞ i=0 (1 -aq i ).
The Heine's summation formula states that for any two real numbers a 1 and a 2 , and for any b with |b/(a 1 a 2 )| < 1, we have

(2.18) 2 ϕ 1 (a 1 , a 2 ; b; q; b/a 1 a 2 ) = (b/a 1 ; q) ∞ (b/a 2 ; q) ∞ (b; q) ∞ (b/a 1 a 2 ; q) ∞ .
We now use Heine's summation formula to conclude the proof: let N be an integer larger than 1; setting a 1 = q -N and a 2 = a yields the q-Chu-Vandermonde formula (2.19) 2 ϕ 1 (q -N , a; b; q; b a q N ) = (b/a; q) N (b; q) N .

Since (q -N ; q) n = 0 as soon as n ≥ N + 1, the series on the left is a finite sum, and we obtain (2.20)

N n=0 (q -N ; q) n (a; q) n (q; q) n (b; q) n b a q N n = (b/a; q) N (b; q) N .
This is an equality between rational fractions in q, which is valid for all q different from 0, 1, and all a = 0 and all b which are not an integer negative power of q. Choosing a = -q and b = -q -N +2 , we have

(2.21) (b/a; q) N = (q -N +1 ; q) N = (1 -q -N +1 ) • • • (1 -q -N +1 q N -1 ) = 0,
and thus the right hand side in (2.20) vanishes, and we get

(2.22) N n=0
(-q; q) n (q -N ; q) n (q; q) n (-q -N +2 ; q) n q n = 0.

We express now p n in terms of q-analogues,

(2.23)

p n = n i=1 1 + q i 1 -q i = n-1 i=0 1 + qq i 1 -qq i = (-q; q) n (q; q) n .
Defining qN-n := (q -N ; q) n (-q -N +2 ; q) n q n and q N -n := qN-n q0 , (2.22) becomes

N n=0 p n qN-n = N n=0 p n q N -n = 0.
Introducing the definition of q N -n and performing similar steps as in (2.23) backwards shows that q n is equal to the value in (2.14), and the proof is complete.

In the steps (a) and (c) of the direct time parallel solver (1.7) based on diagonalization, the condition number of the eigenvector matrix S has a strong influence on the accuracy of the results, and normalizing the eigenvectors from Lemma 2.5 with respect to the 2 norm leads to an asymptotically better condition number, (2.24)

S := V D, D = diag( 1 1 + N -n i=1 |p i | 2
).

We now study the roundoff error when solving the system (B + a 2 I)u = f numerically using the diagonalization

(2.25) S(D + a 2 I)S -1 u = f .
Due to roundoff, we will obtain an approximate solution û, and the difference between the exact solution u and û is classically related to the condition number of the matrix S. We first also give such an estimate, but later provide a more accurate one using the structure of the problem at hand. Lemma 2.6. Let u be the exact solution of (B + a 2 I)u = f , and û be the numerically computed solution using the factored form (2.25), and let u denote the machine precision. Then for any norm,

(2.26) u -û u ≤ cond(B) δB B ≤ (2N + 1)u B -1 |S| |S -1 | D + aI ,
where |S| denotes the matrix containing the elements of S in modulus.

Proof. Using backward error analysis [START_REF] Golub | Matrix computations[END_REF], the computed solution satisfies the perturbed systems

(S + δS 1 )ĝ = f , (D + a 2 I + δD) ŵ = ĝ, (S -1 + δS 2 )û = ŵ,
and since S and S -1 are triangular and D is diagonal we get (see [START_REF] Golub | Matrix computations[END_REF])

|δS 1 | ≤ N u|S| + O(u 2 ), |δS 2 | ≤ N u|S -1 | + O(u 2 ), |δD| ≤ u|D + aI| + O(u 2 ).
Solving numerically (B + a 2 I)u = f using the factored form (2.25) is equivalent to solving exactly (S + δS 1 )(D + a 2 I + δD)(S -1 + δS 2 ) û = f , which is of the form

(B + δB)û = f , |δB| ≤ (2N + 1)u|S| |S -1 | |D + aI| + O(u 2 ).
The relative error then satisfies (see [START_REF] Golub | Matrix computations[END_REF])

(2.27) u -û u ≤ cond(B) δB B ≤ B -1 δB ,
and inserting δB from before gives the result in (2.26).

We next provide an asymptotic estimate for the term |S| |S -1 | on the right in (2.26) in the case of a geometric time partition using the infinity norm. Lemma 2.7 (Asymptotic condition number estimate). For q = 1 + ε, we have

|S| |S -1 | ∞ ∼ 2 2(N -1) (N -1)! ε -(N -1) , (2.28) cond ∞ (S) ∼ 2 (N -1) N N 2 ! N -1 2 ! ε -(N -1) . (2.29) Proof. Note first that |q n | ∼ |p n | ∼ 2 n n! ε n . We next define γ n := 1 + N -n j=1 |p j | 2 and dn := 1
γn , which implies that D = diag( dn ). Then γ n ∼ |p N -n |, and we obtain (2.30)

S ∞ = max n n j=1 |p n-j | γ j ∼ max n n j=1 |p n-j | |p N -j | ∼ max n n j=1 (N -j)! (n -j)! ε 2 N -n ∼ N.
By definition

S -1 = D -1 V -1 = D -1 T (q 1 , • • • , q N -1 ), that is the line n of T (q 1 , • • • , q N -1 ) is multiplied by γ n . Therefore (2.31) S -1 ∞ = max n γ n n-1 j=0 |q j | ∼ max n γ n |q n-1 | ∼ max n γ n |p n-1 | ∼ max n |p N -n ||p n-1 | ∼ 2 ε (N -1) max 1≤n≤N -1 1 (n -1)!(N -n)! .
The last term on the right hand side can be explicitly computed by noting that max 1≤n≤N -1

1 (n -1)!(N -n)! = 1 (N -1)! max 1≤n≤N -1 N -1 n -1 = 1 (N -1)! N -1 N -1 2
, which follows from the Pascal triangle, and thus we get

max 1≤n≤N -1 1 (n -1)!(N -n)! = 1 N 2 ! N -1 2 ! .
Combining (2.30) and (2.31) then yields (2.29). Similarly, we also obtain

|S| |S -1 | ∞ = max i j (|S| |S -1 |) ij = max i i j=1 i k=j |p i-k | |q k-j | ∼ 2 N -1 ε N -1 N k=1 1 (N -k)!(k -1)!
and explicitly summing the last term on the right,

N k=1 1 (N -k)!(k-1)! = 2 N -1
(N -1)! , leads to the desired estimate (2.28). We are now ready to give a precise asymptotic error estimate of the roundoff error that is induced by the diagonalization in the direct time parallel solver (1.7): Theorem 2.8 (Asymptotic roundoff error estimate). Let u be the exact solution of (B + a 2 I)u = f , and û be the computed solution from the direct time parallel solver (1.7) applied to (2.1), and let u denote the machine precision. Then

(2.32) u -û ∞ u ∞ ψ 1 ( aT 2N , N )u ε -(N -1) ,
where

ψ 1 (y, N ) := 2 2(N +1) (N -1)! (1 + 2N (N -1))(1 + y 2 ).
Proof. The form of the matrices involved is very well adapted to estimates in the ∞ norm. We evaluate the various quantities in (2.26) starting with the norm of B -1 , which we can obtain by computing B -1 = (B -1 1 C) 2 explicitly and noting that the infinity norm comes from the last line, which gives

B -1 ∞ = (2k 1 ) 2 1 + 2q(q -1 -1)(q N -1) (q -1) 2 .
Now for q = 1 + ε, we get by expanding

B -1 ∞ ∼ (2k 1 ) 2 (1 + 2N (N -1)).
We now use that D = diag( 4

k 2 i
), and obtain (2.33)

B -1 ∞ D+aI ∞ ∼ (2k 1 ) 2 (4/k 2 1 +a 2 )(1+2N (N -1)) ∼ 16(1+ a 2 k 2 1 4
)(1+2N (N -1)).

Inserting (2.33) into (2.26) gives u -û ∞ u ∞ 16(1 + a 2 k 2 1 4 )(1 + 2N (N -1))u |S| |S -1 | ∞ .
Replacing the last term on the right hand side using (2.28) leads to

(2.34) u -û ∞ u ∞ 2 2(N +1) (N -1)! (1 + 2N (N -1)) 1 + a 2 k 2 1 4 u ε -(N -1) .
Approximating k 1 by k as ε goes to zero finally gives

u -û ∞ u ∞ 2 2(N +1) (N -1)! (1 + 2N (N -1)) 1 + aT 2N 2 u ε -(N -1) ,
which proves the result. Remark 2.9. The more usual round off estimates use the condition number of the matrix S in (2.26) instead of the norm of |S| |S -1 |. Then using (2.29), the estimate (2.32) would become

(2.35) u -û ∞ u ∞ ψ 2 ( aT 2N , N )u ε -(N -1) ,
where ψ 2 (y,

N ) := N (1+2N (N -1)) 2 N +3 N 2 ! N -1 2 !
(1 + y 2 ), which we will see is comparable to our first estimate. We next present a sharper estimate using the special structure of S and S -1 . Lemma 2.10. For any f , for any diagonal matrix ∆, let u = S∆S -1 f , and û be the computed value of u. Then

(2.36) u -u ∞ 2 2(N -1) (N -1)! ε -(N -1) u ∆ ∞ f ∞ .
Proof. Standard truncation error estimates as in Theorem 2.8 show that the approximate value of SDS -1 f gives an error which can be bounded by

|S∆S -1 f -S∆S -1 f | ≤ f ∞ ∆ ∞ u      1 1 + |p 1 | + |q 1 | . . . |p N -1 | + |p N -2 |(|q 1 | + 1) + • • • + |q N -1 | + |q N -2 ||p 1 | + • • • + |q 1 | + 1      +O(u 2 ).
Considering the leading term in ε, suppose that ε and N are such that the largest coefficient p N -1 satisfies |p N -1 |u 1, to obtain

S∆S -1 f -S∆S -1 f ∞ f ∞ ∆ ∞ u N -1 n=0 |p j q N -1-j |.
We can estimate the sum using an asymptotic expansion,

N -1 n=0 |p j q N -1-j | ∼ N -1 n=0 2 n n! ε -n 2 N -1-n (N -1 -n)! ε -(N -1-n) ∼ 2 2(N -1) (N -1)! ε -(N -1) ,
and we obtain (2.37)

S∆S -1 f -S∆S -1 f ∞ f ∞ ∆ ∞ 2 2(N -1) (N -1)! ε -(N -1) u,
which concludes the proof.

We can now prove the following sharper estimate for the roundoff error: Theorem 2.11 (Sharper asymptotic roundoff error estimate). For any f , let u = S(D + aI) -1 S -1 f , and û be the computed value of u. Then

(2.38) u -û ∞ U 0 ψ 3 ( aT 2N , N )u ε -(N -1) ,
where

ψ 3 (y, N ) := 2 2N -1 2 N (N -1)! 1 y 2 +1 . Proof. We have (D + a 2 I) -1 ∞ = 1 a 2 + 4 k 2 N
, and replacing ∆ by (D + a 2 I) -1 in (2.36), we obtain

(2.39) u -u ∞ 2 2(N -1) (N -1)! ε -(N -1) u 1 a 2 + 4 k 2 N f ∞ .
We need now an estimate of f ∞ , with (2.40)

f = C -1 (B 1 C -1 F 0 + G 0 )e 1 , F 0 = 1 k 1 f + g 2 , G 0 = 1 k 1 g - a 2 2 f.
We first explicitly compute the inverse of C and obtain the lower triangular matrix C -1 ij = 2(-1) i-j , i ≥ j. We can then evaluate the terms in (2.40) to get

(C -1 e 1 ) i = 2(-1) i-1 =⇒ (C -1 B 1 C -1 e 1 ) i = 4(-1) i-1   1 k 1 + 2 i j=2 1 k j   ,
where the sum equals zero for i = 1. Using that the time steps are geometric, k j = q j-1 k 1 , we define

s i := 1 + 2 i-1 j=1 1 q j for i ≥ 2, s 1 = 1,
and obtain (C -1 B 1 C -1 e 1 ) i = 4 k1 (-1) i-1 s i , which we introduce into (2.40) to get

f i = (-1) i-1 ( 4 k 1 F 0 s i + 2G 0 ) = (-1) i-1 4 k 1 s i ( 1 k 1 f + g 2 ) + 2( 1 k 1 g - a 2 2 f ) .
Taking the modulus and rearranging terms leads to

|f i | = 4 k 2 1 f s i - a 2 k 2 1 4 + 2 k 1 g(s i + 1) .
Defining y 1 := k1a 2 , we obtain

|f i | = 4 k 2 1 f (s i -y 2 1 ) + 2a k 1 g a (s i + 1) = 4 k 2 1 f (s i -y 2 1 ) + y 1 g a (s i + 1) .
For sufficiently small k 1 , y 1 is smaller than 1, and since s i > 1 and the sequence s i is increasing, we obtain

f ∞ ≤ 4 k 2 1 (|f |s N + |g| a (s N + 1)) ≤ 4 k 2 1 (s N + 1)(|f | + |g| a ) ≤ 4 √ 2 k 2 1 (s N + 1) U 0 .
We use now that q is close to 1 to estimate asymptotically

s N + 1 = 2(1 + N -1 j=1 1 q j ) = 2 N -1 j=0 1 q j = 2 1 -1 q N 1 -1 q 2N,
and therefore

f ∞ 8N √ 2 k 2 1 U 0 .
This allows us to also estimate f ∞ in (2.39) asymptotically,

u -u ∞ 2 2(N -1) (N -1)! ε -(N -1) u 1 a 2 + 4 k 2 N 8N √ 2 k 2 1 U 0 .
We now combine the last two terms asymptotically using that

4 k 2 1 a 2 + 4 k 2 N ∼ 1 1 + y 2 , y := aT 2N ,
and obtain (2.38). We show in Figure 2.2 graphically that ψ 1 ( aT 2N , N ) and ψ 2 ( aT 2N , N ) are comparable, while ψ 3 ( aT 2N , N ) is substantially smaller, and thus gives a much sharper estimate.

Optimization of the algorithm parameters.

To start, we perform a numerical experiment for the scalar model problem (2.1) with a = 1 and initial conditions u(0) = 1, u(0) = 0. We show in Figure 2.3 first as a reference the measured discretization error on an equidistant time partition with N points (called 'error 1', and independent of ε). Next we plot the error due to the geometric time partition, which increases when ε is increasing (called 'error 2'), and our theoretical estimate φε 2 from Theorem 2.1, which is rather sharp, since it coincides in the figure with the numerically measured values, except for large ε. We finally also plot the roundoff error we measure due to the diagonalization procedure, which decreases when ε is increasing (called 'error 3'), and our theoretical bound ψ 3 ε -(N -1) u from Theorem 2.11. From these plots, we can see that an optimized choice of ε would balance these two errors, and for the case on the left with T = 5 and N = 10 (10 processors) it would be ε * ≈ 4.5e -2. In this case, the additional errors due to the time parallelization would remain much smaller than the actual discretization error one would have on an equidistant time grid. For the case on the right in Figure 2.3 with T = 10 and N = 20 (20 processors), the best choice would be ε * ≈ 1e -1, and now the additional errors due to the parallelization would be of the same order as the actual discretization error one would have on an equidistant time grid.

To determine a formula for the best ε * , we can bound the actual error of the parallel time integrator by adding and subtracting u N (T 1 ) and also u N (T q ), with q = 1 + ε, and using the triangle inequality, we obtain at time T an error estimate between the exact solution at time T = t N and the last value ûN of the approximate solution computed by diagonalization, namely

(2.41) |u(t N ) -ûN | U 0 ≤ |u(t N ) -u N (T 1 )| U 0 ≤ error 1 + |u N (T 1 ) -u N (T q )| U 0 ≤ error 2 + |u N (T q ) -ûN | U 0 ≤ error 3
.

The first term error 1 is the truncation error of the sequential method using equal time steps, see also Figure 2.3. The second term error 2 is due to the geometric time partition and was estimated asymptotically in Theorem 2.1 to be φε 2 . The last term error 3 can be estimated using Theorem 2.11. Because the second term is decreasing in ε and the last term is growing in ε, see also Figure 2.3, we equilibrate them asymptotically and obtain Theorem 2.12 (Optimized geometric time partition). Suppose the time steps are geometric, k n = q n-1 k 1 , and q = 1 + ε with ε small. Let be machine precision. Fix a, T and N . For ε = ε * (aT, N ) with

(2.42) ε * (aT, N ) = 3 2 2N (N 2 -1)(N -1)! 1 + y 2 y 3 u 1 N +1
, with y = aT 2N , the error due to time parallelization is asymptotically comparable to the one produced by the geometric time partition. Proof. Equilibration of the error produced by the geometric mesh and the diago- 1) , and thus

nalization means imposing φ(y, N )ε 2 = ψ 3 (y, N ) u ε -(N -
ε * = ψ3 φ (y, N ) u 1 N +1 ,
where ψ 3 is defined in Theorem 2.11, and φ is defined in Theorem 2.1. Introducing these quantities and simplifying leads to (2.42). We see in Figure 2.3 that the theoretically predicted optimized ε * marked by a rhombus is a good estimate of the numerical best parameter. We next show in Figure 2.4 on the left the optimized value ε * (aT, N ) from Theorem 2.12 as a function of the two arguments aT and N . Choosing ε = ε * (aT, N ), the ratio between the additional errors due to parallelization to the truncation error of the fixed time step method (error 2/error 1) is shown in Figure 2.4 on the right, together with the relative truncation error for a fixed time step method in red (error 1). We see that one can use approximately 20 processors (N = 20) for a reasonably large range of aT by only increasing the error with the same order one would have had using a sequential integrator on a fixed time mesh, and when using 10 processors, the additional errors due to time parallelization are negligible (about 10% of the errors on an equidistant time grid). When using 40 processors however (N = 40), a five fold error has to be expected, due to the parallelization, compared to the sequential time stepping on a fixed time mesh, and this when choosing the best possible value ε = ε * (aT, N ). This shows that one can not parallelize with the diagonalization technique in time using an arbitrarily large number of processors.

Remark 2.13. To obtain an estimate for the best ε to chose in the wave equation case, we recall the Fourier transform in space with Fourier variable ξ which corresponds to our parameter a. We can thus apply our results from the ODE analysis with a 2 = |ξ| 2 , where ξ is the dominant frequency in the solution we are trying to calculate. If we have no information about the dominant frequency in the solution, we can also use the upper bound ξ 2 max := π 2 /h 2 in 1D (h the mesh size in space), and |ξ max |2 := π 2 /h 2 1 + π 2 /h 2 2 (h 1 and h 2 the two mesh sizes in space in 2D), since φε * (aT, N ) 2 is growing in aT , see Figure 2.4 on the right 2 . Note that the growth seems to slow down for large values of aT , so that one can even estimate the value of ε * (aT, N ) outside the plotted range of aT .

3. Numerical Experiments. We now test our theoretical results on the wave equation (1.1) in one and two spatial dimensions, and also on an industrial test case using the equations of elasticity, where we use time windowing in order to integrate over a larger number of time steps than authorized by the diagonalization technique.

3.1. The wave equation in one dimension. We solve the wave equation in one dimension on [0, 1], with homogeneous Dirichlet boundary conditions and as initial conditions u 0 (x) = sin(πx), u 1 (x) = 0. The spatial mesh is h = 1 10 , and we use N = 10 time steps on the time interval (0, 1). We first show in Figure 3.1 three numerical solutions obtained for three different geometric time partitions using the diagonalization approach. We clearly see on the left that the ε chosen is too small, the solution becomes inaccurate and presents strong oscillations due to the roundoff error introduced by the diagonalization. On the right ε is too big, and the solution is less accurate for the second part of the time interval. In the middle, we chose the optimal ε * and obtain an accurate solution, comparable to the one computed sequentially on a fixed step size partition. In Figure 3.2 we show the corresponding measured discretization and diagonalization errors compared to our theoretical bounds. We see again on the left that for T = 1 using N = 10 ( 10 time steps corresponding to computing with 10 processors), using the optimized choice of ε leads to a very small error increase compared to the already existing truncation error, as in the ODE case in Figure 2.3. For N = 20 however (corresponding to using 20 processors) on a twice as long time interval T = 2, the errors due to parallelization are now comparable to the existing truncation error on an equally spaced time grid, and as in the ODE case, one can not use more than about 20 processors for the time parallelization of the wave equation without introducing substantial additional errors due to the parallelization method based on diagonalization.

3.2. The wave equation in two dimensions. We now solve the wave equation on the unit square, with homogeneous Dirichlet boundary conditions and as initial conditions u 0 (x, y) = sin(πx) sin(πy), u 1 (x, y) = 0, using h = 1 10 in space. In Figure 3.3 we show again the corresponding measured discretization and diagonalization errors compared to our theoretical bounds, and we see that the results are very similar to the one dimensional case.

We also implemented this algorithm on a computer cluster adding to the initial algorithm a time windowing process in order to compare costs and accuracy for a larger number of time steps, and to truly evaluate parallel performance. Time windowing means that the initial global time interval (0, T ) is first partitioned into so called time windows of length T , and then the diagonalization technique is applied sequentially to the time windows, starting with the first one, then continuing with the second one, and so on. We denote by NWin the number of time windows, where each time window consists of N time steps to apply the diagonalization technique in each time window.

The algorithm now performs 5 steps for each time window, namely

(1)

F = ( 1 k1 u 0 h + 1 2 u0 h , • • • , 0) G = ( 1 k1 u0 h + 1 2 ∆ h u 0 h , 0, • • • , 0) (a) g = (S 1 ⊗ I x )F + (S 2 ⊗ I x )G, (b) ( 1 k 2 n -∆ h )w n = g n , 1 ≤ n ≤ N, (c) u = (S ⊗ I x )w, (2) u = (B -1 1 C ⊗ I x )∆ h u + (B -1 1 ⊗ I x )G,
where

S 1 = S -1 C -1
Y = (A ⊗ I x )X ⇐⇒                    Y 1 = A(1, 1)X 1 Y 2 = A(2, 1)X 1 + A(2, 2)X 2
. . .

Y i = A(i, 1)X 1 + A(i, 2)X 2 + . . . + A(i, i)X i . . . Y N = A(N, 1)X 1 + A(N, 2)X 2 + . . . + A(N, N )X N .
In our implementation, we overlap computation and communication when working in a given time window: once the data from step (1) in (3.1) is available, to perform the products of the form (3.2) in step (a), processor 1 sends its data X 1 to processors 2 to N before computing Y 1 . Processor 1 can then start with its task for step (b). Similarly, processor 2 sends its data X 2 to processors 3 to N , then computes A(2, 2)X 2 and as soon as it receives X 1 also A(2, 1)X 1 to form the sum Y 2 and then starts its step (b). Processors 3 to N proceed similarly. Steps (c) and ( 2) which contain also products of the form (3.2) are performed in the same way overlapping computation and communication. Once the last processor N has completed its step (2), it then becomes processor 1 for the next time window to save communication. Note also that the processors never need to store the entire space time solution, they only need to compute locally the sums in (3.2) using storage of the size of the vectors X i to compute Y i .

In order to measure the we used in this implementation a 1 finite element discretization for step (b), and solved the linear system in space using the conjugate gradient method (CG) with and without preconditioner, and also the direct solver PARDISO [START_REF]Intel MKL PARDISO[END_REF]. Computations were run on a set of 16 X5670 [Westmere-EP] @ 2.93 GHz processors with a refined mesh, h = to T = 2 using NWin time windows of length T = T /NWin. The error is measured in the norm defined in the theoretical analysis, see Theorem 2.1. We show in Table 3.1 the optimized parameter ε * from our analysis given by formula (2.42), and the error due to the time parallelization by diagonalization of the algorithm, which was the same in all experiments, because we used the very small tolerance 10 -16 for the iterative solvers. We can see in the last column that for 16 and more processors, the error due to parallelization will deteriorate, as predicted by our analysis.

In Table 3.2, we show the results we obtain using unpreconditioned CG. The first row corresponds to the sequential scheme with a fixed time step k = T 128 . We then use NWin time windows depending on N to get a result over the entire time interval (0, T ). The column NIt gives the range of the number of iterations CG needed to converge. The column Time RHS corresponds to steps (1) and (a) in algorithm (3.2), the column Time CG to step (b), and the column Time U, Up to steps (c) and (2). The efficiency Eff := T ime1 N ×T ime N , where T ime 1 is the sequential computing time on one processor, and T ime N is the parallel computing time using diagonalization with N processors. Two efficiencies are given in the last two columns: Solv. Eff concerns step (b) and Tot. Eff concerns the whole time sequence. We see from column NIt that as local problems become less well conditioned, the number of iterations increases: we need between 15-33 iterations using one processor, and between 5-149 when using 16 processors, which also affects load balancing. In Table 3.3, we use CG preconditioned with an incomplete LU factorization (ILU0). This improves convergence, but now the cost of an iteration is about twice the cost compared to the non preconditioned case; nevertheless the overall computation time is faster. In Table 3.4, we used a direct solver based on the efficient PARDISO (mkl) library [START_REF]Intel MKL PARDISO[END_REF] that renumbers mesh nodes so that the factorized symmetric matrix is as sparse as possible. We see that PARDISO provides very good local efficiency, and is faster than the iterative methods we used in this 2D setting. This good local efficiency however leads to a little less total efficiency than the variant using CG.

These results show that time parallelization by diagonalization is very efficient as long as the local problem is sufficiently computation intensive so as not to be overtaken by communications cost necessary to compute the right hand side (RHS) and the final solution (U) and its first time derivative (Up).

3.3. Application to a large 3D industrial problem. We now apply this technique to an industrial problem: we want to compute the response of a carbon/epoxy laminated composite panel subjected to an impact-like loading. This class of material is modeled in its linear elastic domain with a transverse isotropic Hooke law. When the fiber direction is aligned with the e 1 direction, and using the Voigt notation3 , the Hooke law is given by

(3.3)         ε 11 ε 22 ε 33 2ε 23 2ε 13 2ε 12         =          1 E L -ν LT E L -ν LT E L 0 0 0 -ν LT E L 1 E T -ν T T E T 0 0 0 -ν LT E L -ν T T E T 1 E T 0 0 0 0 0 0 2 1+ν T T E T 0 0 0 0 0 0 1 G LT 0 0 0 0 0 0 1 G LT                  σ 11 σ 22 σ 33 σ 23 σ 13 σ 12         .
Here, E L is the longitudinal elastic modulus (in the carbon fiber direction), E T is the transverse elastic modulus, G LT is the elastic shear modulus, and ν LT and ν T T are the so-called Poisson ratios, and we simulate the elasticity equations, ρü = div(σ) + f , ε ij = 1 2 ∂u i ∂x j + ∂u j ∂x i .

The selected material is a T700GC/M21 carbon/expoxy system used in the aeronautical industry. The corresponding elastic moduli are given in Table 3.5. The geometry of the problem is taken from the experimental work described in [START_REF] Trousset | Prediction of low-velocity impact damage on composite plates[END_REF]. The tested preprocessing and the MPI communications are shown separately, and we see that the LU linear solves are fully parallelized in time with this diagonalization technique. Figure 3.5 indeed shows that the corresponding times scale nicely up to N = 8 processors, and one obtains substantial speedup using this method. The overhead due to the different algebraic operations (mainly factorization and tensorial operations) remains about the same from N = 2 to N = 16 but it becomes more significant in the last case (steps (a), (c) and (2) in (3.1)). From a practical point of view, the 650s gain from a sequential to an N = 8 parallel solution seems to be narrow. However, this kind of simulation is typically used to optimize the mechanical behavior according to the stacking sequence. Depending on the optimizer, this process may not be fully parallelizable. When tens or hundreds of runs are required, the gain of using time parallelization then becomes significant. Figure 3.6 shows the displacement of the central node on the back face of the plate for the sequential solution and the parallel solution for N = 16. No deviation is noticeable between the two responses, showing that no significant error has been introduced by the parallelization with these parameters.

This type of time parallel algorithm would be well suited to a hybrid MPI/OpenMP approach: MPI for space domain partitioning which can by far exceed shared memory systems and OpenMP for the time parallelisation (less than 20 threads). In terms of memory, there could be a gain (space variables would not need to be replicated) but probably not so in terms of CPU time which is already quite good as shown in the histogram in Figure 3.5 (see MPI gathering and MPI scattering).

4. Conclusions. We presented and analyzed the time parallelization method by diagonalization for the wave equation. We derived an optimal choice for the geometric stretching of the time grid, balancing carefully truncation and roundoff error. We tested the method numerically for the model wave equation problems in different spatial dimensions, and also on an industrial test case using the elasticity equations. These results show that our theoretical parameters predict well the best choice for the geometric time grid stretching, and substantial speedup is possible when solving wave propagation problems using this technique, and this also in an industrial setting.

. 15 ) 2 n

 152 Proof. The matrix B = (C -1 B 1 ) 2 is diagonalizable if and only if all time steps k n are different. The eigenvalues are then 4 k , and the eigenvectors are those of C -1 B 1 .

Fig. 2 . 2 .Fig. 2 . 3 .

 2223 Fig. 2.2. Comparison of the logarithm of the functions ψ j , j = 1, 2, 3

Fig. 2 . 4 .Fig. 3 . 1 .

 2431 Fig.2.4. Left: optimized choice ε * (aT, N ) from Theorem 2.12. Right: ratio of the additional errors due to parallelization to the truncation error of the fixed step method with this choice of ε * (aT, N ) (error 2/error 1), and relative truncation error for a fixed time step method in red (error 1)

Fig. 3 . 2 .Fig. 3 . 3 .

 3233 Fig. 3.2. Discretization and parallelization errors in 1d, together with our theoretical bounds for the PDE. Left: T = 1, N = 10. Right: T = 2, N = 20

  B 1 and S 2 = S -1 C -1 . Compared to the three basic steps (a), (b) and (c) of the initial algorithm (1.7), the additional step (1) serves to set up the proper initial values of the current time window, and the additional step (2) has to be computed in order to initialize F and G for the next time window.All the time matrices are lower triangular matrices and all Kronecker products Y = (A ⊗ I x )X where X = (X 1 , X 2 , ..., X N ) and Y = (Y 1 , Y 2 , ..., Y N ) can be written as(3.2)

Fig. 3 . 5 .Fig. 3 . 6 .

 3536 Fig. 3.5. Computing times for the industrial elasticity problem.

Table 3 . 1

 31 Optimal ε * and common error to all solvers N NWin NIt Time RHS Time CG Time U, Up Solv. Eff Tot. Eff

	1	128	15-33	0.001 s	0.780 s	0.057 s
	2	64	15-38	0.025 s	0.460 s	0.041 s	84.78 %	79.69 %
	4	32	14-47	0.037 s	0.269 s	0.059 s	72.49 %	57.40 %
	8	16	10-63	0.045 s	0.160 s	0.139 s	60.94 %	30.45 %
	16	8	5-149	0.064 s	0.101 s	0.180 s	48.27 %	15.18 %

Table 3 .

 3 2 2D wave equation using CG without preconditioning, initialization time 0.0007 s

N

  NWin NIt Time RHS Time-CG-ILU Time U, Up Solv. Eff Tot.

								Eff
	1	128 11-11	0.001 s	0.665 s	0.062 s
	2	64	11-11	0.026 s	0.383 s	0.041 s	86.81 % 80.89 %
	4	32	10-11	0.038 s	0.166 s	0.071 s	100.15 % 58.39 %
	8	16	7-15	0.051 s	0.088 s	0.119 s	94.46 % 33.09 %
	16	8	2-33	0.059 s	0.074 s	0.129 s	56.17 % 17.37 %
					Table 3.3
		2D wave equation using CG and ILU0 preconditioner, initialization time 0.0039 s, including
	ILU0 factorization
	N NWin NIt Time RHS Time Direct Sol. Time U, Up Solv. Eff Tot. Eff
	1	128	1	0.001 s	0.256 s	0.033 s
	2	64	1	0.024 s	0.134 s	0.034 s	95.52 % 64.36 %
	4	32	1	0.038 s	0.066 s	0.067 s	96.97 % 37.76 %
	8	16	1	0.047 s	0.035 s	0.130 s	91.43 % 17.10 %
	16	8	1	0.061 s	0.016 s	0.146 s	99.99 %	8.13 %

Table 3 .

 3 4 2D wave equation using PARDISO, initialization time 0.102 s including the factorization with PARDISO

Nick Higham, after a presentation of the first author, and later also an anonymous referee, suggested considering a random partition. We tested this and did not find any advantage over the geometric partition in our numerical experiments. The geometric partition allows us to provide a complete analysis for the best choice of the parameters in the method, while for the analysis of random time steps it is not clear how one would have to proceed.

In fact we need only spectral information of the spatial operator to determine the relevant range of a, and thus can also handle variable coefficients, see the example in Subsection

3.3.

Voigt notation uses a matrix representation for

th order tensors. It also exploits the symmetries of the 2 nd order strain and stress tensors (hence the 6 components).

plate is made of 12 plies with a symmetric stacking sequence. The mesh, illustrated in Figure 3.4, consists of 152607 degrees of freedom. Each ply is modeled with one element in the thickness direction. The stacking sequence (i.e. the angle between the fiber orientations and the global e x direction for each ply) is taken as [ 0 / 45 / 90 / -45 / 0 / 0 / 0 / 0 / -45 / 90 / 45 / 0] which implies that the material is heterogeneous at the structure scale. Dirichlet boundary conditions (represented in blue) are applied on areas representative of the experiment. A homogeneous pressure field is applied at the center of the plate (red circle) and is representative of a small impact. The time is discretized with 2000 time steps over the T :=10ms simulation range.

The finite element code Zset [START_REF]Z-set, material and structure analysis suite[END_REF] is used to handle the specific finite element features: mesh generation, material orientation, computation of mass and stiffness matrices, boundary conditions, etc. Zset comes with a Python interface which permits to handle the global operators (stiffness and mass matrices) and vectors in a very convenient way through numpy arrays. All the algebraic operations described in equations (1.2) to (1.7) are thus performed with the numpy and scipy Python packages. Parallelism is achieved with the mpi4py package, the Python interface to MPI. The script is launched on a 2.4GHz, Xeon E5-2680 processor with 24 cores, but we use only up to N = 16 for our experiments. To reach the 2000 time steps, we thus have to use 2000/N time windows. There is an evident gain here in using a direct solver, since once the N sparse operators (see equation (1.7.b)) have been factorized during the first block, their decomposition can be reused for the solutions in the 2000/N -1 remaining time windows. In this implementation, we use the MUMPS solver [START_REF]MUMPS: a MUltifrontal Massively Parallel sparse direct Solver[END_REF] and its LU decomposition capability. Because of the multiple time windows, we also have to properly account for the non-zero initial conditions at the beginning of each block. These operations require additional algebraic operations that we refer to as windows initializations (steps (a), (c) and ( 2) in (3.1)).

We show in Figure 3.5 the measured CPU times when parallelizing with different values of N . The times for the LU factorizations, the LU linear solves by forward and backward substitution, the tensorial operations, the windows initializations, the