Luis T Antelo

Stéphanie Passot

Fernanda Fonseca

Cristian Ioan

Trelea

Antonio A Alonso

Towards optimal operation conditions of freeze-drying processes via a multi-level approach

Keywords: Freeze Drying, Modelling, Dynamic Optimization, Control Profiles, Real Time Optimization

Freeze-drying (lyophilization) offers an attractive dehydration method for valuable food and biological products, as it is capable of preserving product quality and biological activity while extending their shelf life. However, despite these benefits in terms of product quality, freeze-drying is also a notoriously energy-intensive and time-consuming process. This requires an expensive operation to construct an efficient optimal decision-making tool able to drive the operation through the most effective paths, that minimize time and maximize product quality. Here we propose an integrated approach to operational design and control of the freeze-drying process that combines dynamic modelling with efficient optimized off-line and on-line control. The required mass and energy balance equations still contain inherent nonlinearity, even in their lumped parameter version. This results in a set of complex dynamic, computationally-costly optimization problems solved by selected global stochastic optimization algorithms. Real-time disturbances and model uncertainties are addressed via the proposed hierarchical multi-level approach, allowing recalculation of the required control strategies. The framework developed has been revealed as a useful tool to systematically define off-line/on-line optimal operation policies 0 for many food and biological processing units.

Introduction

Freeze-drying is a stabilization process typically used in food and biological industries. It provides the highest quality dried product by freezing the material and heating it to force the frozen water to transition directly from solid to gas phase (sublimation), producing an interconnected porous structure that can be rehydrated very quickly. Sublimation is an endothermic phase transition that, for pure substances, occurs at temperatures and pressures below the triple point. Compared to classical thermal drying methods, freezedrying better preserves the biological activity of thermosensitive components as well as the organoleptic and nutritional properties of the material. Lyophilized products have a longer shelf-life and are more convenient for transport and storage. The main disadvantage of this process stems from slow drying rates due to working conditions that require low pressure and temperature, making freeze-drying a very expensive process in terms of time, energy, and therefore economic costs. This limits the use in the food industry, being only viable for dehydration of high added-value products such as pharmaceuticals and living cells.

The end-to-end freeze-drying process comprises three main steps: freezing, primary drying, and secondary drying. The principal function of the freezing process is to separate the solvent from the solutes. In an aqueous system, the water will form ice crystals while the solutes will be confined to the interstitial region between ice crystals. Once the product has reached a completely frozen state, the pressure in the freeze-dryer is reduced, and heat is applied to the product to initiate sublimation of the ice crystals. The sublimation process gradually creates a moving interface between frozen zone and dried zone. During primary drying, which is the longest part of the process, product temperature needs to be kept below the collapse temperature. Primary drying is complete once all ice crystals have been removed from the product, and the volume occupied by the resulting cake is equivalent to that of the frozen matrix. On completion of sublimation, there will be some water adsorbed onto the cake surface. This moisture can constitute up to 5-10% of the dried product depending on the temperature and nature of the cake components. In many cases, these moisture values may prove too high, and the final product may not have the desired stability. The desired stability is obtained during secondary drying by desorbing the moisture from the cake without reducing its interstitial volume. Removal of the unfrozen water is usually accomplished by increasing the product temperature while reducing the pressure. Typical figures for final moisture level are less than 3%.

As cited in [START_REF] Trelea | An interactive tool for the optimization of freeze-drying cycles based on quality criteria[END_REF], various process variables affect the efficiency of the freeze-drying process. Product temperature is a determinant factor for both productivity and product quality. Indeed, sublimation is generally faster at higher temperatures, but drying at excessive temperatures results in a loss of the pore structure obtained by freezing, which in [START_REF] Pikal | The collapse temperature in freeze-drying: Dependence on measurement methodology and rate of water removal from the glassy state[END_REF] is defined as the collapse phenomenon and therefore in rejection of the batch. The collapse temperature is usually about 2 °C higher than the glass transition temperature of the maximally freeze-concentrated phase (T Glass ) or equal to the eutectic temperature if solutes are crystallized in the frozen solution. As presented by [START_REF] Alves | Advances in multi-purpose drying operations with phase and state transitions Drying[END_REF], the appropriate shelf temperature and chamber pressure conditions are frequently established empirically in a trial-and-error experimental way. In this framework, the simplest approach to determine optimal operation policies consists in considering constant values of these variables, which are generally set on the basis of results obtained in previous runs carried out with the same product. As a result, this approach usually leads to non-optimal conditions for the freezing, sublimation and desorption steps [START_REF] Liapis | Research and development needs and opportunities in freeze-drying[END_REF]. Significant advantages can be obtained if shelf temperature and chamber pressure are varied during the operation [START_REF] Pisano | In-line optimization and control of an industrial freeze-drying process for pharmaceuticals[END_REF]. Therefore, several methodologies proposed in the literature are based on employing process simulation to determine the operating conditions over the shelf temperature and chamber pressure that minimize drying time while satisfying the glass transition constraint [START_REF] Rene | Vacuum freeze-drying of a liquid in a vial: determination of heat and mass-transfer coefficients and optimisation of operating pressure[END_REF]- [START_REF] Gan | Lyophilization in vials on trays: Effects of tray side on lyophilisation performance[END_REF]. Proper definition of the mathematical model provides a better understanding of process dynamics and their influence on cycle time and product temperature history [START_REF] Liapis | Freeze-drying of pharmaceutical crystalline and amorphous solutes in vials -Dynamic multidimensional models of the primary and secondary drying stages and qualitative features of the moving interface[END_REF]- [START_REF] Hottot | Freeze-drying of pharmaceutical proteins in vials: modeling of freezing and sublimation steps[END_REF], which is crucial for obtaining truly optimal process operation that delivers high quality product at minimal cost [START_REF] Sadikoglu | Freeze-drying of pharmaceutical products: Research and development needs[END_REF].

A systematic procedure to optimally determine time-varying operating profiles is one of the scopes of recent research on freeze-drying [START_REF] Sadikoglu | Optimal control of the primary drying of solutions in vials using variational calculus[END_REF] - [START_REF] Fissore | In-line control of a freeze-drying process in vials[END_REF], [START_REF] Pisano | In-line optimization and control of an industrial freeze-drying process for pharmaceuticals[END_REF]. In this context, the present work aims to go a step further by considering dynamic optimization to define a multi-level integrated approach for optimal control of food and biotechnological processing units.

The tool developed distinguishes between an a priori (off-line) optimization level based on the defined model and a real time (on-line) optimisation level. The off-line optimization level will set up nominal optimal operation policies for the process, while the on-line level takes into account unpredicted disturbances and unavoidable model uncertainties that might affect the process in order to recalculate valid the operational profiles. These profiles are introduced into a low level regulatory layer (actuators) which forces plant states to evolve accordingly.

The aim of the approach presented in this paper is to systematically determine optimal operation policies for freeze-drying processes, while attempting to overcome the disadvantages detected in previous approaches. Several studies have reported one stage off-line optimization strategies, producing a fixed process operation recipe [START_REF] Alves | Advances in multi-purpose drying operations with phase and state transitions Drying[END_REF], [START_REF] Gan | Lyophilization in vials on trays: Effects of tray side on lyophilisation performance[END_REF], [START_REF] Sadikoglu | Optimal control of the primary and secondary drying stages of bulk solution freeze drying in trays[END_REF], [START_REF] Millman | An analysis of the lyophilization process using a sorption-sublimation model and various operational policies[END_REF]. Their main drawback is the sub-optimal character of the resulting operation policies, either because of unexpected disturbances, model parameter uncertainty or equipment variability [START_REF] Pisano | In-line optimization and control of an industrial freeze-drying process for pharmaceuticals[END_REF]. Sub-optimal policies can violate important process constraints as pointed out in [START_REF] Sadikoglu | Optimal control of the primary and secondary drying stages of bulk solution freeze drying in trays[END_REF], [START_REF] Sadikoglu | Optimal control of the secondary drying stage of solutions in vials using variational calculus[END_REF]. For example, they may lead to product collapse (local melting) which in turn can put at risk the entire batch. In contrast, the proposed multi-level optimisation approach can efficiently prevent those violations by re-computing optimal operating policies in real time based on actual process measurements.

The off-line level used here to analyze different processing scenarios defined as a function of product moisture content, is described in Section 2. The main results obtained for the different control configurations are summarized in Section 3. Finally, Section 4 presents the principle of the hierarchical approach based on implementing the precomputed policies in a real-time optimization (RTO) framework designed to react to possible disturbances and minimize the adverse effects of final product variability.

The off-line level: The freeze-drying optimal control problem

The aim of this level is to set up optimal operation policies in time (t) over considered controlled variables (u) by solving an optimization problem of the form:

min u J s.t. f(ż, z, p, u, t) = 0 h(z, p, u, t) = 0 g(z, p, u, t) ≥ 0 u l ≤ u ≤ u u (1)
where J is a given performance index, z ∈ R α is the state vector; ż is the time derivative of z; p ∈ R np is a set of possible process parameters; u ∈ R nu is the control vector; f is the set of differential and algebraic equations describing the system's dynamics and h and g are possible path, point or end point constraints (equality and/or inequality), which represent further special conditions governing process operation. The search space for the control variables u is delimited through upper and lower bounds (u l and u u ).

For the freeze-drying case considered in this work (schematized in Figure 1), the objective is to define process operation policies which minimize cycle duration, i.e. J = t cycle , by varying the defined control variables u= [T shelf , P t chamb , T cond , T chamb ] in order to attain a given final average moisture content in the product (C ave ) while meeting product stability conditions, which in this case means that product temperature (T ) has to remain below the glass transition temperature (T Glass ) at any point in the entire freeze-drying cycle. Mathematically:

min u J = t cycle s.t. f(ż, z, p, u, t) = 0 C Ave (t f inal ) = 0.03 kg water /kg T ≤ T Glass -50 ≤ T shelf ≤ 30(°C); 0 ≤ P t chamb ≤ 50(Pa); -90 ≤ T cond ≤ -50; 5 ≤ T chamb ≤ 50 (2)
where the target average moisture content of the product (C Ave ) is considered as an endpoint constraint while the quality requirements (T ≤ T Glass ) are managed as path constraints. Present work integrates the one-dimensional dynamic model (f ) of heat and mass transfer with associated product quality indicators developed by [START_REF] Trelea | An interactive tool for the optimization of freeze-drying cycles based on quality criteria[END_REF] for the freezedrying process of PS product (4% of polyvinylpyrrolidone (PVP), 1% of sucrose, and 10 mM Tris-HCl). The complete differential-algebraic model equations (f ) considered are reproduced in full in Appendix 1.

It should be stressed that in practice T cond and T chamb are not usually considered as manipulated variables, but the aim of our work is to explore all control alternatives/strategies available that could potentially help improve freeze-drier design and performance. Therefore, different control scenarios based on the definition of vector u will be analyzed:

1. One control variable: u = [T shelf ].

Two control variables

: u = [T shelf , P t chamb ].

Four control variables

: u = [T shelf , P t chamb , T cond , T chamb ].
Furthermore, note that the lower bound for the controlled variable P t chamb (0 Pa) is not physically possible, but this limiting case has been examined in order to see whether lowering this bound would result in a significant improvement in terms of freeze-drying cycle duration. Later in this work, a more realistic bound of 10 Pa will be introduced for this variable.

These different process configurations will be translated in particular versions of the previously defined dynamic optimization problem (2). To solve them, a Control Vector Parametrisation approach ( [START_REF] Vassiliadis | Computational solution of dynamic optimizations problems with general differantial-algebraic constraints[END_REF], [START_REF] Vassiliadis | Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints[END_REF] -Figure 2) is used to transform the original dynamic optimization problem of infinite dimension into a non-linear optimization problem (NLP) of finite dimension. For this purpose, the considered time interval is divided into ρ constant time intervals and the controls are approximated in each interval generally by loworder polynomials, as depicted in Figure 2. This new problem can be solved by employing different optimization algorithms considering that process dynamics need to be integrated on each internal iteration in order to evaluate both the objective function and constraints (if any). In this framework, considering high ρ values would entail substantial computational costs without significantly improving the final solution of the performance index J [START_REF] García | Computing optimal operating policies for the food industry[END_REF]. On top of these important solving requirements associated to higher-dimension problems (these problems used to be NP-Hard, meaning that the computation time required to solve them increases exponentially with their dimension), higher discretization levels lead most optimization approaches to suboptimal solutions with very noisy control profiles [START_REF] Balsa-Canto | Dynamic optimization of distributed parameter systems using second-order directional derivatives[END_REF], [START_REF] Chachuat | Global methods for dynamic optimization and mixed-integer dynamic optimization[END_REF]. In the other hand, overly low ρ values could result in significantly sub-optimal profiles as it would become impossible to obtain sufficient curvature for the control profiles. Therefore, in this work, the number of time intervals considered for approximating the control variables using the CVP approach is ρ = 10 and 20. Larger values of this parameter could be chosen, approximating the control variable to a real continuous profile. However, and as mentioned, the computation time associated to solving an optimization problem of higher dimension could be unaffordable. Moreover, there are many cases to show that global solutions obtained at selected levels of discretization have a very low margin of error (for the case of problems with known solutions) when compared to more refined solutions [START_REF] García | Computing optimal operating policies for the food industry[END_REF], [START_REF] Balsa-Canto | Dynamic optimization of distributed parameter systems using second-order directional derivatives[END_REF]. In this context, the influence of the control discretization (ρ) in the optimum freeze-drying cycle duration will be analyzed later in this work.

From what has been discussed so far, it can be concluded that the complexity of the NLP problems defined for each operational scenario is intimately related to:

1. The highly non linear nature of the freeze drying process and its associated mathematical model which involves a high number of parameters and exhibits a rich dynamics.

2. The existence of constraints over the process states and controls.

3. The number of selected controls.

what in many instances may lead to local optima. Strongly zigzagging control profiles were found for this problem when using non-robust optimizers (FMINCOM, DIRECT, etc.) what can be considered as an evidence of sub-optimal solutions (note that usually, global optimal profiles are expected to be smooth).

As a consequence, optimal control problems are computationally costly as they require:

• A numerical integration scheme to solve the set of ordinary differential equations within each objective function evaluation.

• To apply methods to approximate the selected controls by polynomial functions (as the CVP approach considered in this work) in order to transform the original dynamic optimization problem of infinite dimension into a non-linear optimisation problem (NLP) of finite dimension.

In this off-line framework , and to overcome the above mentioned drawbacks, three global stochastic optimisation algorithms have been considered for solving the low-dimension resulting NLP problems (derived from 2):

• DE: Differential Evolution [START_REF] Storn | Differential evolution -a simple and efficient heuristic for global optimization over continuous spaces[END_REF]. This is a heuristic algorithm for the global optimization of nonlinear and (possibly) non-differentiable continuous functions. It is a population-based method which, starting with a randomly generated population, computes new candidate solutions by calculating differences among population members. It handles stochastic variables by means of a direct search method which outperforms other popular global optimization algorithms, and it is widely used by the evolutionary computation community.

• G-CMA-ES: Covariance Matrix Adaptation Evolutionary Strategy [START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF]. This is an evolutionary algorithm that makes use of the covariance matrix in a similar way to the inverse Hessian matrix in a quasi-Newton method, and it is particularly interesting for solving ill-conditioned and non-separable problems. This method was ranked in the first place in the Conference of Evolutionary Computation 2005 (CEC'05).

• eSS-SSm: Enhanced Scatter Search [START_REF] Egea | Scatter search for chemical and bio-process optimization[END_REF]. This method was recently developed for solving nonlinear dynamic optimization problems, outperforming other stateof-the-art methods. It is based on the use and combination of a reference set of good solutions, maintaining an appropriate level of diversity.

The rationale for these selections is that stochastic methods usually provide excellent solutions in affordable computation times. The state-of-the-art in deterministic methods for dynamic optimization problems is currently unsatisfactory, as in practice, methods exploiting the model structure can only be used for very small problems, otherwise the computational times to achieve the global solutions become unaffordable.

Regarding deterministic methods, several of them can be employed (e.g., a multistart procedure using a SQP -Sequential Quadratic Programming, method, a branch and bound method or similar). However, if the problem is multimodal, multistart SQP may fail to locate the global optimum (as shown in next Section), whereas the use of a branch and bound method is again restricted to small problems. Research comparing stochastic and deterministic methods has shown that the stochastich algorithms are more efficient in terms of computation time [START_REF] Egea | Scatter search for chemical and bio-process optimization[END_REF], [START_REF] Banga | Dynamic optimization of bioprocesses: Efficient and robust numerical strategies[END_REF] and [START_REF] Egea | Dynamic optimization of nonlinear processes with an enhanced scatter search method[END_REF].

In order to select which of these algorithms leads to the best results, it is necessary to perform an efficiency analysis. For this purpose, convergence curves are constructed showing the evolution of the best value obtained by each solver as a function of CPU time. With these representations, both the robustness (the capability of the solver to attain consistently good final solutions) and the efficiency (the speed to converge to the final solution) can be evaluated, allowing the user to select the appropriate algorithm for solving a given optimization problem, as presented later in this work.

Results on dynamic optimization of freeze-drying

This Section presents the solutions of the resulting dynamic optimization problem (2) for the different control scenarios considered. These solutions have been obtained by using the selected optimization algorithms.

One control variable: T shelf profile

As a first step, the same scenario discussed in previous works on freeze-drying optimization [START_REF] Trelea | An interactive tool for the optimization of freeze-drying cycles based on quality criteria[END_REF] is considered by solving problem (2) when only one control variable is taken into account (u= [T shelf ]). In this case, P t chamb = 10 Pa, T cond = -65 °C and, T chamb = 25 °C are constant along the equidistant ρ=20 intervals in which the considered time interval (t = 40 h) has been split.

In order to demonstrate the complexity of the defined optimal control problem as well as the inability of local deterministic algorithms to attain the global optimum, a multistart strategy with a SQP method (FMINCON from MATLAB © [START_REF]The Mathworks. Optimization Toolbox 4 User guide[END_REF]) has been performed by considering 30 random vectors for T shelf with a control discretization of ρ = 20.

As shown in the histogram of achieved solutions given in Figure 3 for the multistart strategy, the problem is non-convex (multiple different solutions are obtained depending on the T shelf profile considered) and, therefore, stochastic global optimization algorithms are required to properly solve it.

It needs to be pointed out that when using stochastic optimization methods like the ones considered in this study, it is necessary to perform several runs to increase the probability of finding the global optimum. Depending on the computational cost of the problem, the literature usually cites numbers between 10 and 30. Therefore, and as a consequence of the high complexity and computational cost associated to the defined optimal control problem, a number of 10 independent runs for each algorithm tested (performed on a Intel Core2 Quad 2.40 GHz CPU with a fixed maximum computation time of 8 hours) was chosen. We will consequently obtain an average on the freeze-drying cycle time.

The results achieved for the optimum cycle duration when using eSS-SSm (hereafter denoted SSm), DE and G-C-MAES are presented in Table 1. As shown on this Table, all average values for the selected objective function (freeze-drying cycle duration) are fairly similar, with the one obtained by SSm being slightly better (J Ave = 25.067 h) than the ones reached by DE and C-MAES (J Ave = 25.076 h and 25.180 h, respectively). In a real-world process, these differences would not be significant when compared to other operations such as freeze-drier loading, sterilization, freezing, etc. Note that Table 1 also shows that C-MAES exhibits a higher (but not very significant) variability (measured through the standard deviation) on the final value of J obtained for each run. Either way, all these optimizations obtain better solutions than the best one achieved by the multistart strategy with FMINCON (J = 27.02 h).

As previously cited in Section 2, deterministic global algorithms will fail to reach the global optimum due to the highly non-linear multimodal nature of problem (2). To illustrate this fact, this simplest optimal control problem has been solved with DIRECT (DIviding RECTangles, [START_REF] Jones | DIRECT global optimization algorithm[END_REF]), which is an algorithm based on a modification of the Lipschitzian optimization scheme. It operates by systematically dividing the optimization domain in hyper-rectangles and evaluating the objective function in their centres. The optimum freeze-drying cycle time obtained by DIRECT is J = 28.02 h, which is close to 8% higher than the best value obtained by SSm. As a consequence, this type of algorithms is no longer considered in this work.

As stated above, in order to properly select a reference solver, it is necessary to perform an efficiency analysis based on the constructed convergence curves (Figure 4). First of all, we can check that all the runs performed for each solver converge to substantially similar final J valued. In Figure 5, solvers efficiencies are compared by considering the convergence curve corresponding to the best run obtained by each (Run 6 for SSm and DE and 2 for C-MAES). The so-called best run is the one that provides the best solution in terms of objective function within the maximum computation time allowed. It obviously provides a high-quality solution even if it is impossible to demonstrate that it is certainly the global optimum. Stochastic methods such as those employed in this work do not provide proofs for convergence, but they usually provide excellent solutions in affordable computation times.

For this scenario, the optimum control profile (corresponding to T shelf ) obtained by the different algorithms are shown in Figure 6. It can be concluded from this figure that the achieved solutions are almost identical, leading to fairly similar cycle durations.

Comparing the presented results against those reported in [1] (optimum cycle duration -J best = 24.4 h), it can be seen that this value is lower than the best solution (J best = 25.067 h) obtained by solving the considered dynamic optimization problem [START_REF] Pikal | The collapse temperature in freeze-drying: Dependence on measurement methodology and rate of water removal from the glassy state[END_REF]. This difference could be due to the fact that in [START_REF] Trelea | An interactive tool for the optimization of freeze-drying cycles based on quality criteria[END_REF], the product top temperature is higher than the corresponding glass transition temperature at the beginning of the freezedrying cycle, thus violating the product quality requirement, as shown in Figure 7. In this setting, when the optimum T shelf profile obtained in [START_REF] Trelea | An interactive tool for the optimization of freeze-drying cycles based on quality criteria[END_REF] (Figure 7) is compared with those obtained in this work via dynamic optimization with a path constraint (Figure 6), the main difference between them is the T shelf value in the second control discretization point. Such difference guarantees, for the proposed optimal profile, that T top ≤ T top Glass throughout the cycle (since T shelf is lower), preventing the system from collapsing but also slightly increasing the freeze-drying process duration.

Two control variables: T shelf and P t chamber profiles

This part of the work analyzes the influence of changes on chamber pressure over the freeze-drying cycle and its duration, this time considering two control variables to define the NLP problem presented in ( 2 As with one-control variable case, several runs were carried out for the three algorithms selected to solve the NLP problem presented in (2) considering the new defined vector u. The results obtained using SSm, DE and C-MAES are presented in Table 2. In this particular case, and as explained in Section 2, two levels of control discretization have been considered (ρ = 10 , 20) in order to determine the influence of this parameter on the optimum freeze-drying cycle time. The reason for this selection has been already put forward in Section 2: high ρ values would entail substantial computational costs without significantly improving J, while overly low discretization levels could result in significantly sub-optimal profiles as it would be impossible to obtain sufficient curvature for the control profiles. In order to maintain a reasonable compromise between both extremes, affordable discretization levels of 10 and 20 can be selected.

In order to experimentally prove these statements, a set of tests have been carried out by solving [START_REF] Pikal | The collapse temperature in freeze-drying: Dependence on measurement methodology and rate of water removal from the glassy state[END_REF] with SSm for different values of ρ = [START_REF] Pisano | In-line optimization and control of an industrial freeze-drying process for pharmaceuticals[END_REF][START_REF] Liapis | Freeze-drying of pharmaceutical crystalline and amorphous solutes in vials -Dynamic multidimensional models of the primary and secondary drying stages and qualitative features of the moving interface[END_REF][START_REF] Velardi | Development of simplified models for the freeze-drying process and investigation of the optimal operating conditions[END_REF][START_REF] Egea | Scatter search for chemical and bio-process optimization[END_REF]40,50]. For each discretization level, ten runs were performed with a fixed maximum computation time of 24 hours. Such limit was selected to be of the order of the optimal freeze-drying cycle time. The reason is that in order for the algorithm to be applicable in real time, what implies re-calculation (and implementation) of optimal policies in the event of disturbances affecting the process, computation cannot take longer than the duration of the real process itself. As a result, we obtain the J Ave as a function of ρ, as depicted in Figure 8. This Figure shows how the average freeze-drying cycle time quickly decreases when increasing ρ up to values between 20 and 30, at which point J Ave reaches a minimum before starting to slightly increase for higher values of control discretization. Therefore, it can be concluded that the computational effort tied to high ρ values does not significantly improves the final solution. In fact, medium levels of discretization (ρ = [START_REF] Liapis | Freeze-drying of pharmaceutical crystalline and amorphous solutes in vials -Dynamic multidimensional models of the primary and secondary drying stages and qualitative features of the moving interface[END_REF][START_REF] Velardi | Development of simplified models for the freeze-drying process and investigation of the optimal operating conditions[END_REF]) can be considered representative for optimal control and real time optimization purposes.

From the results presented in Table 2, it can be concluded that:

1. A higher control discretization (ρ = 20) leads to better J values, i.e to lower lyophilization cycle times, for all the solvers considered but entails a strong increase in computational effort to achieve the optimum. This can be expected since with a higher discretization level, the parametrized control profile becomes a better approximation of the theoretical continuous profile.

2. The final values of J obtained for ρ = 20 clearly improve (by a 5-6%) the values obtained by all solvers under the one-control-variable case. This means that the freeze-drying cycle duration is cut by almost an hour and a half (from J Ave = 25.067 h. to J Ave = 23.506 h.). Note that these values correspond to the best solutions obtained by SSm. As a result, the influence of changes in P t chamber profiles over the cycle duration is demonstrated.

3. The average J values obtained by SSm are better than those obtained by DE and C-MAES for both levels of discretization considered.

Focusing back on the efficiency of the solvers, Figure 9 shows that the different runs performed by the solvers considered only converge to fairly similar final J values for the SSm case, while C-MAES and especially DE lead to higher dispersion in final solutions, revealing a low robustness of these solvers when solving (2) for this two-control variable scenario. This fact is corroborated by the high values for the standard deviation of the mean value presented in Table 2 for CMAES and DE (σ=[0.227,0.4649], respectively) compared to SSm (σ=0.060).

In Figure 10, the efficiencies of the solvers are compared by considering the convergence curve corresponding to the best run (leading to the minimum J) obtained by each. This aspect will prove to be very important when the on-line optimization level is defined. In this context, the solvers employed should be "sufficiently fast" to react against possible disturbances and to minimize the adverse effect of final product variability but in lower computational times, as detailed later. In this framework, it can be seen from Figure 10 that CMAES and SSM are revealed as good options for real time optimization purposes.

Finally, the optimum control profiles (corresponding to T shelf and P t chamb ) obtained are shown in Figure 11. It is concluded that, for this case, the profiles are quite different, with evident deviations for the total pressure chamber case. This situation might indicate the existence of a number of local minima that are very similar in terms of the objective function value but somewhat different in terms of decision variables. However, the global shape of the control profiles is similar, indicating that the algorithms converged to physically meaningful profiles.

Two control variables: Ramps in T shelf and steps in P t chamber profiles

Practical considerations for the operation of freeze-dryers on real food and biotech processing plants make it necessary to introduce a new operational issue regarding the employment of step changes (piecewise constant intervals for control discretization) over P t chamber . This is due to the fact that this variable is difficult to modify/manipulate a real scale, making it necessary to implement profiles based on step changes.

From this concrete operational scenario, and as a result of what has been stated in the previous cases, SSm has been considered as the most suitable algorithm for solving the resulting dynamic optimization problem. Once again, two levels of control discretization have been considered (ρ = 10 , 20) in order to analyze the influence of this parameter on the cycle duration. When comparing the values presented in Table 2 for step profiles in the chamber pressure (Table 2, fourth column) to the values obtained when considering ramps in both manipulated variables (Table 2, third column), it can be concluded that the average of the cycle duration (J) is slightly increased when using steps in the P t chamber profile, while reducing the standard deviation of the mean value of J and, therefore, variability in the solutions attained in different runs.

Regarding the influence of ρ, Figure 12 compares the convergence curves for the best runs performed for each control discretization level, showing, as in Table 2, that the absolute minimum for the freeze-drying cycle duration is obtained for a higher level of discretization (ρ = 20), but at the price of a higher (around one order of magnitude) computational effort. Once again, a compromise between robustness and efficiency is needed when defining the real time optimization (RTO) level, since the J value obtained for ρ = 10 could be good enough and the time needed to achieve it much lower (around one order of magnitude) than for the ρ = 20 case.

The optimum control profiles (corresponding to T shelf and P t chamb ) obtained with SSm are depicted in Figures 11, 13 and 14. These Figures make possible to check that the temperature profiles for the control discretization levels considered (ρ) are smoother than those obtained when ramps are used for both control variables, as shown in Figure 11a. This makes easier to translate/apply these profiles into real-world operation of the freezedrying process.

Finally, it must be pointed out that real equipment limitations make it impossible to achieve in-chamber pressures under 7-8 Pa to generate total vacuum. As a consequence, the lower bound for the P t chamb control variable has been modified from 0 to 10 Pa (P t chambu ≥ 10) with a discretization level ρ = 10. The average value of the cycle duration is almost the same as the one obtained for the case of considering unconstrained P t chamb (J Ave = 24.681 h versus J Ave = 24.671 h, value extracted from Table 2, ρ = 10).

These optimum profiles obtained when solving the constrained P t chamb two-controlvariables DO problem are presented in Figure 14 and compared with the previous nolower-bound case solution. As shown in this Figure, the temperature profile for this constrained scenario is almost the same as that obtained for the unconstrained pressure case, while the pressure profile varied as expected when compared to the unconstrained case due to the new set of the lower bound for P t chamb ,. The penalization on the final freeze-drying cycle time is almost negligible.

Four control variables: T shelf and P t chamber profiles, T cond and T chamb

The last configuration analyzed for the freeze-drying process considers four control variables: T shelf and P t chamber profiles, T cond and T chamb . Note that only T shelf and P t chamb are parametrized along the homogeneous partition of the time interval [0, 40] into ρ subintervals.

As mentioned earlier, this operational scenario can not be practically implemented on most of the current real freeze-driers, but the aim of this work is to investigate and analyze possible advantages that could shape future equipment design.

The results obtained when solving the complete NLP problem presented in Equation (2) using SSm, which has been revealed as the more appropriate solver (in terms of efficiency-robustness), conclude that the average value achieved for the freeze-drying cycle duration is J Ave = 23.332 h. This means that the effect of changes on T cond and T chamb only slightly modify optimal operation time, since a reduction of only 0.17 h (a 0.074%) is attained when compared to the two-control-variables scenario (for the higher control discretization level considered -ρ = 20). In short, the gain is minor given the effort of modifying freeze-dryer design, thus reinforcing current real-scale freeze-drying equipment configurations.

Figure 15 presents the optimum control profiles (corresponding to T shelf and P t chamb ) for this four-control variables case. In addition to these profiles, optimum values for nondiscretized decision variables have been obtained. Optimal values for T cond and T chamb are T condopt = -87.219 °Cand T chambopt = 26.347 °C, respectively.

To summarize this Section, a comparison among the optimum profiles obtained for all cases considered (one, two and four-control variables) is presented in Figure 16. It can be concluded that the trend in the optimum T shelf profiles is the same to the trend reported in [START_REF] Trelea | An interactive tool for the optimization of freeze-drying cycles based on quality criteria[END_REF], since during primary drying, shelf temperature has to be decreased because the selfcooling effect due to ice sublimation decreases when mass transfer resistance through the dry layer increases. Subsequently, and during secondary drying, the product moisture content decreases and the glass transition temperature increases throughout in the product. This leads to the possibility of significantly increasing the shelf temperature. Looking at the chamber pressure profiles, the trend along the considered time interval is similar for both considered cases (only ramp variations have been compared).

4 Re-computing control profiles on disturbed freeze-drying processes: a first step to define a RTO level Many food and biotech industries use time-constant control profiles on their processes, but allowing the control profiles to change in time can improve process performance. As described earlier in this paper, dynamic optimization (DO) uses mathematical models of the processes to compute optimal control profiles that maximize/minimize a predefined performance index such as productivity or other economic indexes.

Usually, these control profiles are computed off-line and implemented in the plant. However, when an unpredicted disturbance (for instance, changes in the external temperature, input/output fluid pump failure, etc.) enters the plant, the control profiles computed using the DO scheme are no longer optimal. A possible framework to over-come this limitation is real-time optimization (RTO), in which the control profiles are re-computed a number of times throughout the process.

The general scheme of RTO is presented in Figure 17. Sensors collect information on the state variables of the real plant. This information is fed into the optimization software which, via the simulation environment (mathematical model), computes new optimal control profiles. Finally, these profiles are introduced into a low-level regulatory layer (actuators) which will force the states of the plant to evolve accordingly.

In Section 2, we saw that the objective of dynamic optimization of freeze-drying is to minimize the cycle duration by varying the defined control variables in order to attain a given average moisture content in the final product while satisfying product stability conditions related to glass transition temperature.

However normal operation of freeze-drying process could still be affected by different disturbances or uncertainties (both in the model parameters as well as in the measured or manipulated variables). To overcome their negative effects on the process dynamics, an RTO level can be defined to repeatedly re-compute the control profiles obtained through dynamic optimization (and which are no longer optimal). With the aim of defining a proper and reliable RTO environment, a set of simulation tests have been carried out to represent typical operation scenarios under the existence of disturbances. The objective is to determine how fast the selected optimization solvers (the core of the RTO level) can lead the system to a new optimal operation point by recalculating control profiles. Mathematically, the resulting problem is like (4), but now including the vector of possible disturbances affecting the process during the entire time interval considered (d), and using off-line optimum profiles (u * ) as starting point. This work considered three step disturbances affecting the system dynamics in order to define the freeze-drying operation case studies proposed:

1. A +10% step on T cond (From -65 to -58.5 °C) 2. A +10% step on T chamb (From 25 to 27.5 °C) 3. A -10% step on T chamb (From 25 to 22.5 °C)

From the analysis of the results achieved for dynamic optimization of freeze-drying cycles, global solvers emerged as the best alternatives to re-computing control profiles. Once again, three selected global NLP solvers have been employed to solve this new optimization problem: SSm, DE and CM AES.

In addition, and with the aim of possibly improving RTO level efficiency, we also considered hybrid strategies, methods that combine global stochastic and local optimization algorithms by taking advantage of both the robustness of stochastic solvers and the efficiency of local methods when started in the optimum neighbourhood. Here, three hybrid strategies were: a) SSm+FMINCON (local gradient-based method included in the Matlab Optimization toolbox - [START_REF]The Mathworks. Optimization Toolbox 4 User guide[END_REF]); b) SSm+MISQP (an extention of the well-known SQP method to mixed-integer programmes - [START_REF] Exler | A trust region sqp algorithm for mixed-integer nonlinear programming[END_REF] ; and c) SSm+IPOPT (interior point optimizer for large-scale nonlinear optimization - [START_REF] Wächter | On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming[END_REF]).

To avoid overloading the reader with an array of numerical and graphical results, Table 3 simply presents the best, worst and mean values of J (freeze-drying cycle duration) achieved for each solver for the operation scenarios considered under the presence of a step disturbance on T cond (+10%). All solvers attain good solutions for the new optimum J (in the neighbourhood of 25 h) in the interval 0.15-0.40 h. (≈ 10-25 min.), the best value of the objective function being obtained by SSm (J best eSS = 24.665 h.) while DE posted the worst performance of all algorithms considered. To verify these results, Figure 18 presents the convergence curves corresponding to the best run (in terms of final value of J) for each solver/hybrid strategy considered. These results indicate that the proposed algorithms, which constitute the basis of the RTO level, are able to re-compute the disturbance-rejecting control profiles that alter optimum process operation in a time that can be considered minimal given the global freeze-drying cycle duration (≈ 1 -2% of the total FD cycle time).

If this on-line level is not considered, and the optimal profiles obtained off-line (as presented in Figure 14) are used to try to reject the considered disturbance in condenser temperature, then we get a violation of the glass transition temperature constraint. This is due to the fact that the disturbance knocks process operation away from the optimal point determined off-line, making it an unfeasible operating scenario with an associated cycle time that increases up to J = 31.14 h. Clearly, re-computing control profiles via RTO is a justified measure for overcoming this drawback.

It is also stated in Figure 18 that, for this specific case study, employing hybrid strategies does not actually increases the efficiency of the global solvers since no better solutions are found in shorter times compared to global solvers such as CMAES or SSm. It can be also concluded that SSm, whether as a global solver or as part of an hybrid strategy, consistently achieves high-quality solutions (best J value) and thus present the best compromise between efficiency and robustness.

Finally, Figure 19 presents the comparison between the off-line optimum control profiles and the re-computed control profiles when a step on T cond is considered. It shows that the profiles for the shelf temperature are almost the same for the process with and without disturbance, while pressure profiles reflect significant changes designed to reject the disturbance. Similar results are obtained when other considered disturbances affect the process, as stated in Figure 20.

Conclusions

This work presents a multi-level integrated approach that defines operation conditions for minimizing freeze-drying cycle time while preserving product quality (final water content) by solving a dynamic optimization problem. The high dimensionality and nonlinearity inherent to the model equations describing the full freeze-drying dynamics result in a set of computationally-involved control problems which justify this two-level hierarchical approach.

The off-line level is responsible for analyzing different processing scenarios. In this framework, several scenarios have been considered based on the decision variables selected. This kind of assessment will be formally stated as a dynamic optimization (DO) problem to provide the given average moisture content in the product while satisfying both process dynamics and product stability and quality conditions. By solving the different DO problems with the global stochastic solvers considered (SSm, DE and CMAES), it was clearly stated that by considering two control variables (T shelf and P t chamb ) makes it possible to achieve a reduction of up to 1.5 h in the freeze-drying cycle duration compared to the one-control-variable scenario with T shelf . The non-practical four-control-variables DO problem helped us to show that the major effort required to potentially modify and improve the design of today's freeze-driers by also manipulating T cond and T chamb will be not translated into significant process time reduction (no more than a 0.1% gained).

Furthermore, it was also verified that higher levels of control discretization (ρ) leads to better values of cycle times values, achieving a minimum of J ave in the interval between ρ= 20 and ρ=30 before starting to slightly increase at higher values of control discretization. Therefore, it is concluded that the computational effort associated to high values of ρ does not significantly improves the final solution. In fact, medium levels of discretization (ρ = [START_REF] Liapis | Freeze-drying of pharmaceutical crystalline and amorphous solutes in vials -Dynamic multidimensional models of the primary and secondary drying stages and qualitative features of the moving interface[END_REF][START_REF] Velardi | Development of simplified models for the freeze-drying process and investigation of the optimal operating conditions[END_REF]) can be considered representative for optimal control and real time optimization purposes.

We saw that the on-line level consists of implementing the pre-computed policies on a real time optimization (RTO) framework designed to react to possible disturbances and to minimize the adverse effects of end-product variability by re-computing the off-line control profiles. In this framework, and in order to define the fundamentals of this level, several tests were performed to re-optimize the freeze-drying operation under the presence of a step disturbance on T cond (+10%). It was shown that all solvers achieve good solutions (SSm giving the best performance) for the new optimum J (in the neighbourhood of 25 h) in the 0.15-0.40 h time interval (≈ 10-25 min.) following the application of the disturbance. At this point, it can be concluded that a RTO level based on the considered optimization algorithms will be able to re-compute the control profiles (which reject the disturbance) in a short time when compared to total freeze-drying cycle duration (≈ 1 -2% of the total freeze drying cycle time).

The authors are currently engaged in research aimed at completing the results presented here by: a) analyzing model uncertainty and its effects on the solution of the optimal control problem; b) defining and implementing a robust RTO level to re-compute these control profiles at given time intervals under disturbances and model uncertainty, and c) implementing off-line computed profiles and the proposed multi-level approach at pilot/real scale.

Appendix 1: Heat and mass transfer model equations for the freeze-drying process

This research considers the model described in [START_REF] Trelea | An interactive tool for the optimization of freeze-drying cycles based on quality criteria[END_REF]. It consists of a system of four differential equations (for the front position and the residual moisture contents at the product bottom, front, and top) and one algebraic equation for the front temperature. Note that the product and the dryer state are described by relevant state variables at six key points, as shown in Figure 1: shelf, product bottom, sublimation front, product top, freeze-drying chamber, and condenser. The heat and mass transfers are assumed to take place between these points.

Primary drying (I)

Secondary drying (II)

Control variables T shelf , T chamb , T cond , P t chamb State variables Z f ront , T f ront , C bottom , C f ront , C top

Mass transfer

Mass transfer coefficients

K F T o = k F T o A Zt-Z f ront K T o Ch = k T o Ch A 1 K F Co = 1 K F T o + 1 K T o Ch + 1 K Ch Co

Heat transfer

Heat transfer coefficients

H S B = h S B (P t chamb , H 2 O)A H S B = h S B (P t chamb , N 2 )A H B F = h B F A Z f ront H F T o = h F T o ( P F w +P T o w 2 , H 2 O) A Z t -Z f ront H F T o = h F T o ( P F w +P T o w 2 , N 2 ) A Z t H T o Ch = h T o Ch A 1 H S F = 1 H S B + 1 H B F 1 H S Ch = 1 H S B + 1 H F T o + 1 H T o Ch 1 H F Ch = 1 H F T o + 1 H T o Ch
Heat fluxes

Q S F = H S F (T f ront -T shelf ) Q S Ch = H S Ch (T chamb -T shelf ) Q F Ch = H F Ch (T chamb -T f ront ) Temperatures T bottom = T shelf + 1 H S B Q S F T bottom = T shelf + 1 H S B Q S Ch T f ront = T bottom T top = T f ront + 1 H F T o Q F Ch

Product quality

Glass transition

T B Glass = T Glass (C B ) T F Glass = T Glass (C F ) T T o Glass = T Glass (C T o )
Sorption isotherms

C B Equ = C Equ (1) C B Equ = C Equ P B w P Sat w (T bottom ) C F Equ = C Equ (1) C F Equ = C Equ P F w P Sat w (T f ront ) C T o Equ = C Equ P T o w P Sat w (Ttop)

State evolution equations

Front position

dZ f ront dt = F F C D Ice A dZ f ront dt = 0 Front temperature 0 = Q S F -Q F Ch -L Subl F F C 0 = T bottom -T f ront

Residual water content

dC bottom dt = 1 τ Des C B Equ -C bottom dC f ront dt = 1 τ Des C F Equ -C f ront dC top dt = 1 τ Des C T o Equ -C top

Additional outputs

Neutral gas press P Ch n = P t chamb -P Ch w Average water content

C Ave = Z F (C bottom +C f ront )/2+(Z t -Z F r )(C f ront +C top )/2 Zt C Ave = C bottom +C top 2
Product mass in a vial

M V ial = M Dry + C Ave 1-C Ave M Dry + Z F Zt M Ice Ini M V ial = M Dry + C Ave 1-C Ave M Dry Functions Sorption isotherm C Equ (a w ) = min C M ax , max C M in , q Sorp 1 (a w -q Sorp 2 )
Transition from primary to secondary drying

f (I) (Z) = sin π 2 Z Z T rans if Z < Z T rans 1 if Z ≥ Z T rans Glass transition T Glass (C) = max T F r Glass , ( 1
N V ial 213 R V ial 7.12e -3 [m] M Dry 0.05e -3 [kg] M Ice Ini 0.95e -3 [kg] Z t 6.82e -3 [m] Z T rans 1.70e -3 [m]
Ice properties

T Ice Glass -135 [°C] D Ice 920 [kg/m 3 ] L Subl 2.83e 6 [J/kg] h BF 2.4 [W/(mK)]
Mass Transfer

k F T o 1.2e -8 [kg/(smP a)] k T oCh 8e -6 [kg/(sm 2 P a)] k ChCo 5e -7 [kg/(sP a)]
Heat transfer between shelf and product bottom

q SB H 2 O 1 q SB N 2 0.625 h SB Co 0.6 [W/(m 2 KP a)] P SB T rans 100 [P a]
Heat transfer between sublimation front and product top

h F T o W alls+Rad 50e -3 [W/(mK)] q F T o H2O 1 q F T o N2 0.813 h F T o Co 327e -6 [W/(mKP a)] P F T o T rans 100 [P a]
Heat transfer between product top and chamber

h T oCh 5 [W/(m 2 K)] Bottom Front Top

  ): u = [T shelf , P t chamb ]. Once again, T cond = -65 °C, T chamb = 25 °C are constant during the defined time interval (t = 40 h).

Figure 1 :

 1 Figure 1: Simplified scheme of a freeze-drier. The controlled variables considered in this work are chamber temperature and pressure (T chamb and P t chamb ) , shelf temperature (T shelf ) and condenser temperature (T cond ).

Figure 2 :

 2 Figure 2: Schematic view of the control vector parametrization applied to a given nonlinear function using zero (steps) and first order (ramps) polynomials.

Figure 3 :

 3 Figure 3: Histogram of solutions attained for the freeze-drying cycle time by the multistart procedure with FMINCON for the one control variable scenario (T shelf ) and ρ=20.

Figure 4 :

 4 Figure 4: Convergence curves obtained with the different solvers considered by solving the dynamic optimization of a freeze-drying cycle. One control variable (T shelf ) and ρ = 20.

Figure 5 :Figure 6 :

 56 Figure 5: Convergence curves for the best run obtained with each solver by solving the dynamic optimization of a freeze-drying cycle. One control variable (T shelf ) and ρ = 20.

Figure 7 :Figure 8 :

 78 Figure 7: Optimum T shelf profile obtained by Trelea et al. [1]

Figure 9 :

 9 Figure 9: Convergence curves obtained with the different solvers considered by solving the dynamic optimization of a freeze-drying cycle. Two control variables: T shelf and P t chamber ; ρ = 20.

Figure 10 :

 10 Figure 10: Convergence curves for the best run obtained with each solver for the dynamic optimization of a freeze-drying cycle. Two control variables: T shelf and P t chamber ; ρ = 20.

Figure 11 :

 11 Figure 11: Optimum profiles obtained for the dynamic optimization of a freeze-drying cycle. Two control variables: T shelf and P t chamber ; ρ = 20.

Figure 12 :

 12 Figure 12: Convergence curves obtained with SSm for the dynamic optimization of a freeze-drying cycle. Two control variables: Ramps in T shelf and Steps in P t chamber

Figure 13 :

 13 Figure 13: Optimum control profiles obtained with SSm for the dynamic optimization of a freeze-drying cycle. Two control variables case with ramps in T shelf and steps in P t chamber and ρ = 20.

  lower bound P t chamber (b) P t chamber

Figure 14 :

 14 Figure 14: Optimum control profiles obtained with SSm for the dynamic optimization of a freeze-drying cycle. Two control variables case with ramps in T shelf and steps in P t chamber and ρ = 10.

Figure 15 :

 15 Figure 15: Optimum profiles obtained with SSm for the dynamic optimization of a freezedrying cycle. Four control variables: T shelf and P t chamber Profiles, T cond and T chamb ; (ρ = 20). The optimum values obtained for T cond and T chamb are T cond opt = -87.219 °Cand T chambopt = 26.347 °C

Figure 16 :

 16 Figure 16: Optimum profiles obtained for the dynamic optimization of a freeze-drying cycle when different number of control variables are considered (ρ = 20).

Figure 17 :Figure 18 :

 1718 Figure 17: General structure of the real time optimization framework

Figure 19 :

 19 Figure 19: Re-computed profiles obtained by using SSm to solve the operation scenario of a freeze-drying cycle with a step disturbance of +10% on T cond . ρ = 10

Figure 20 :

 20 Figure 20: Re-computed profiles obtained by using SSm to solve the operation scenarios of a freeze-drying cycle under different step disturbances. ρ = 10

Table 1 :

 1 Results of solving optimal control problem (2) representing the operational scenario with one control variable: T shelf Profile

	NLP Solver

Table 2 :

 2 Results of solving optimal control problem (2) representing the operational scenario with two control variables: T shelf and P t chamber Profiles

	NLP Solver

Table 3 :

 3 Results for the re-computing of control profiles for the considered freeze-drying process under a step on T cond (+10%)

	NLP Solver

Acknowledgments

The authors received financial support from the 7th EU Framework Programme (CAFE Project -Large Collaborative Project: KBBE-212754), the Spanish Government (MICINN Project AGL2008-05267-C03-01) and Xunta de Galicia (IDECOP 08DPI007402PR). Dr. Luis T. Antelo is grateful for financial support provided under the "Anxeles Alvarino" Programme -Xunta de Galicia.


Alexandre Bellis

On links between horocyclic and geodesic orbits on geometrically infinite surfaces

We study the topological dynamics of the horocycle flow h R on a geometrically infinite hyperbolic surface S. Let u be a non-periodic vector for h R in T 1 S. Suppose that the half-geodesic u(R + ) is almost minimizing and that the injectivity radius along u(R + ) has a finite inferior limit Inj(u(R + )). We prove that the closure of h R u meets the geodesic orbit along un unbounded sequence of points gt n u. Moreover, if Inj(u(R + )) = 0, the whole half-orbit g

Yet, we give a construction where Inj(u(R + )) > 0 and g R + u ⊂ h R u, which also constitutes a counter-example to Proposition 3 of [Led97].

Introduction

Among the curves of constant curvature in H 2 are the geodesics of curvature zero and the horocycles of curvature one. They give rise to two flows with deep relations in the unitary tangent bundle T 1 H 2 : the geodesic flow and the horocycle flow respectively. Consider now a fuchsian group Γ and the quotient surface S := Γ\H 2 . Both of these flows descend to the quotient T 1 S := Γ\(T 1 H 2 ). We denote them by g R and h R respectively.

While orbits of the geodesic flow are known to be diverse, on the opposite, these of the horocycle flow tend to be rigid. It is illustrated by a result of G. Hedlund in [START_REF] Hedlund | Fuchsian groups and transitive horocycles[END_REF] stating that when the surface S is compact, for every u in T 1 S, the orbit h R u is dense in T 1 S. In [START_REF] Eberlein | Horocycle flows on certain surfaces without conjugate points[END_REF], P. Eberlein proves that this result comes from a fundamental link between horocyclic and geodesic orbits : h R u is not dense in the non-wandering set Ω h of the horocycle flow if and only if the projection on S of g R + u, denoted by u(R + ), is almost minimizing (specifically : ∃C > 0, ∀t > 0, d(u(0), u(t)) > t -C).

As a consequence, he obtains that if the surface is geometrically finite (i.e. with a finitely generated fundamental group), if u is in Ω h , then h R u is either dense in Ω h or periodic. This rigidity property was generalised by M. Ratner to Lie groups and unipotent actions in [START_REF] Ratner | Raghunathan's topological conjecture and distributions of unipotent flows[END_REF]. However, it does not extend to geometrically infinite surfaces (i.e. not geometrically finite). Indeed, S is geometrically finite if and only if every horocyclic orbit in Ω h is dense in Ω h or periodic (see [START_REF] Dal | Geodesic and Horocyclic Trajectories[END_REF]).

In this paper, we are interested in the topological dynamics of the horocycle flow on geometrically infinite surfaces, where little is known. Untwisted hyperbolic flutes are the simplest examples of such surfaces and still contain many interesting behaviours (see [START_REF] Haas | Dirichlet points, garnett points, and infinite ends of hyperbolic surfaces. i[END_REF] or [START_REF] Coudène | Horocycles récurrents sur des surfaces de volume infini[END_REF]). More precisely, we investigate the links between the geodesic flow and the horocycle flow. We associate to x in S the real number Inj(x) defined as the maximal radius of a ball centered at x without self-intersection.

The main result of this paper is :

Theorem 1.1. Let Γ be a fuchsian group without elliptic elements and such that the quotient surface S := Γ\H 2 is geometrically infinite. Consider u in the non-wandering set Ω h of h R in T 1 S. Suppose that u(R + ) is almost minimizing and that h R u is not periodic and define Inj(u(R + )) := lim inf t→+∞ Inj(u(t)).

If Inj(u(R + )) < +∞, then there exists a sequence of times t n converging to +∞ such that

g tn u ∈ h R u for all n. Moreover, if Inj(u(R + )) = 0, then g R + u ⊂ h R u.
In particular, when Inj(u(R + )) < +∞ for every u in T 1 S (O. Sarig calls such a surface weakly tame, see [START_REF] Sarig | The horocycle flow and the laplacian on hyperbolic surfaces of infinite genus[END_REF]), then, if h R u is not periodic, there always exists a positive time t such that g t u ∈ h R u. As a corollary, we get : Corollary 1.2. Let Γ be a fuchsian group with no elliptic nor parabolic element such that the quotient surface S := Γ\H 2 is geometrically infinite. If S is weakly tame, then the horocycle flow does not admit any minimal set on T 1 S. This corollary gives an easy way to construct surfaces without a minimal set for the horocycle flow. The first example of such a surface was given by M. Kulikov in [START_REF] Kulikov | The horocycle flow without minimal sets[END_REF]. Later, theorems of non-existence were obtained in [START_REF] Matsumoto | Horocycle flows without minimal sets[END_REF] and [START_REF] Gaye | Sur l'inexistence d'ensembles minimaux pour le flot horocyclique[END_REF].

In the setting of Theorem 1.1, we can ask : [START_REF] Starkov | Fuchsian groups from the dynamical viewpoint[END_REF], A.N. Starkov gives an example of a surface S satisfying the hypotheses of Theorem 1.1 such that 0 < Inj(u(R + )) < +∞ and h R u is not recurrent. In [START_REF] Ledrappier | Horospheres on abelian covers[END_REF], F. Ledrappier generalizes this example to manifolds with bounded geometry (see Proposition 3 of [START_REF] Ledrappier | Horospheres on abelian covers[END_REF]) : let M be a manifold with bounded geometry and u in T 1 M such that u(R + ) is asymptotically almost minimizing. Then, for every t ∈ R, the strong stable leaf

Question : is it possible to have 0 < Inj(u(R + )) < +∞ and g R + u ⊂ h R u, as in the case Inj(u(R + )) = 0 ? Clearly, if h R u is not recurrent (i.e. it does not accumulate on itself), then g R + u ⊂ h R u. This implies in particular that h R u is locally closed (there exists a neighborhood V of u such that V ∩ (h R u -h R u) = ∅) even though it is not closed. In
W ss (g t u) := {v ∈ T 1 M | d(g t+s u, g s v) -→ s→+∞ 0} is locally closed.
When M is a hyperbolic surface S, this statement becomes : suppose S has an injectivity radius everywhere above some positive constant C. Let u be in Ω h such that u(R + ) is almost minimizing. Then h R u is locally closed.

Actually, this proposition is false and I construct in section 4 the following counter-example which also answers the previous question : Theorem 1.3. There exists a geometrically infinite surface S with an injectivity radius everywhere above some positive constant C (i.e. S has bounded geometry), with u in Ω h satisfying :

(i) u(R + ) is almost minimizing. (ii) g R + u ⊂ h R u. (iii) Inj(u(R + )) < +∞.
In particular, h R u is not locally closed.

Acknowledgements. I thank my Ph.D advisor Françoise Dal'Bo for her numerous and precious advices about the redaction of this paper. I also thank Yves Coudène for the fruitful discussions.

Notations and tools

For two points z and z ′ in H 2 and two points ξ and η in ∂H 2 := R ∪ {∞}, we denote by [z, z ′ ] the hyperbolic segment between z and z ′ , by [z, ξ) the halfgeodesic going from z to ξ and by (η, ξ) the geodesic between η and ξ. We denote by g R and h R the geodesic, respectively horocycle flow in the unitary tangent bundle T 1 H 2 . For any u in T 1 H 2 , the function symbol u(t) refers to the projection on H 2 of g t u and u + refers to the extremity in ∂H 2 of the half-geodesic u(R + ).

Consider now two elements u and v in T 1 H 2 . Suppose that u(0) = z and v(0) = z ′ and that u + = v + = ξ. Then there exists t in R such that

g t h R u = h R v. The Busemann cocycle B ξ (z, z ′ ) centered at ξ between z and z ′ is by definition the number B ξ (z, z ′ ) = t. Thus, the set {z ′ | B ξ (z, z ′ ) = 0} is the horocycle centered at ξ passing through z.
Level sets of isometries : The group PSL 2 (R) acts by orientation preserving isometries on (H 2 , d). Let γ ∈ PSL 2 (R) be a hyperbolic isometry. We denote by γ -and γ + the points in ∂H 2 that are its repulsive and attractive extremities respectively. Observe that a point z ∈ H 2 is moved by γ along the hypercycle passing through γ -, z and γ + , that we denote by C γ (z). For a positive integer d, the point γ d z belongs to the portion of C γ (z) between z and γ + .

Let l γ := inf z∈H 2 d(z, γz) be the translation length of γ, realised on its axis (γ -, γ + ). We have : Proposition 2.1 (Section 5 of [START_REF] Parkkonen | On the hyperbolic orbital counting problem in conjugacy classes[END_REF]). For any hyperbolic isometry γ and any

z in H 2 , sinh d(z, γz) 2 = cosh s sinh l γ 2 where s = d(z, (γ -, γ + )).
When γ is parabolic, we denote by C γ (z) the horocycle centered at the unique fixed point γ + = γ -of γ and passing through z.

We have :

Proposition 2.2 (Section 6 of [PP15]
). Consider a parabolic isometry γ and pick any z 0 in H 2 . Denote by l γ (z 0 ) the distance l γ (z 0 ) = d(z 0 , γz 0 ). For any z in H 2 , we have :

sinh d(z, γz) 2 = e s sinh l γ (z 0 ) 2 where s = B γ + (z, C γ (z 0 )).
To prove our theorems, we will translate the dynamics of the horocycle flow on Γ\(T 1 H 2 ) in terms of the action of Γ on H 2 and ∂H 2 using the following proposition.

Proposition 2.3 (Proposition 2.1, Section V of [START_REF] Dal | Geodesic and Horocyclic Trajectories[END_REF]). Take a vector u in T 1 S and a positive real number r. Note ũ a lift of u in T 1 H 2 and suppose that ũ+ is not fixed by any element of Γ except the identity. Then

g r u ∈ h R u -h R u ∃(α n ) N ∈ Γ N s.t. α n ũ+ -→ n→+∞ ũ+ and B ũ+ (α -1 n i, ũ(0)) -→ n→+∞ B ũ+ (i, ũ(r)) .
3 Proof of theorem 1.1

Let us first give a precise definition of the injectivity radius.

Définition 3.1 (Injectivity radius). Let S := Γ\H 2 be a hyperbolic surface. The injectivity radius of S at p is :

Inj(p) := min γ∈Γ γ =id d(p, γ p)
where p is any lift of p to H 2 .

Proof of Theorem 1.1. Consider a lift ũ in T 1 H 2 of u. Up to conjugacy, we can suppose ũ+ = ∞ and ũ(0) = i. Note that with our choice of lift ũ of u, the equivalence of Proposition 2.3 becomes :

(g r u ∈ h R u -h R u) ⇐⇒ ∃(α n ) N ∈ Γ N , α n ∞ → ∞ and B ∞ (α -1 n i, i) → r (1)
The key of this proof is to find elements α n in Γ on which to apply this equivalence.

Step 1 : There is a sequence of points (q n ) n going to ∞ on the half-geodesic [i, ∞) = [ũ(0), ũ+ ) and a sequence of elements γ n in Γ that are all different such that

(i) d(q n , γ n q n ) -→ n→+∞ Inj(u(R + )).
(ii) For every sequence of positive integers

(k n ) N we have γ kn n ∞ -→ n→+∞ ∞.
Proof of Step 1. The hypothesis Inj(u(R + )) < +∞ of Theorem 1.1 gives us a sequence of points q n going to ∞ in ∂H 2 along the half-geodesic [i, ∞) and a sequence of elements

(γ n ) n≥0 in Γ -Id satisfying condition (i) of Step 1. Since neither g R u nor h R u is periodic, these elements γ n are all different.
Consider a subsequence of (γ n ) N , that we still denote by

(γ n ) N , such that lim n→+∞ γ - n = η and lim n→+∞ γ + n = ξ. Suppose first that η = ξ. In this case, for any z in H 2 , the distance s n = d(z, (γ - n , γ + n )
) is bounded from above. Thus, since l γn ≤ d(q n , γ n q n ) for every n, Proposition 2.1 implies that for any z in H 2 , the elements γ n z stay in a compact. This contradicts the discreteness of Γ.

Suppose now that η = ξ = ∞. If the elements γ n are hyperbolic, consider the half-geodesic [p, ξ) starting from a point p on (0, ∞) and orthogonal to (0, ∞) (if η = ξ = 0, consider any point p on (0, ∞)). For n big enough, we have

d(p, (γ - n , γ + n )) < d(q n , (γ - n , γ + n ))
. Thus, Proposition 2.1 implies that d(p, γ n p) < d(q n , γ n q n ). Since the latter is bounded from above, we get again a contradiction with the discreteness of Γ. Finally, if the elements γ n are parabolic, for any z and any z 0 in H 2 , we eventually have B γ + n (z, C γn (z 0 )) < B γ + n (q n , C γn (z 0 )). Thus, Proposition 2.2 implies that d(z, γ n z) < d(q n , γ n q n ), which gives again a contradiction with the discreteness of Γ.

In conclusion,

η = ξ = ∞. (2)
Choose now the following orientation for the elements γ n .

• If γ n is hyperbolic, choose |γ - n | ≤ |γ + n |. • If γ n is parabolic, choose it such that |γ n ∞| > |γ + n |.
This choice of orientation combined with (2) concludes Step 1.

Step 2 : Fix ǫ > 0 and an interval I of R + of length Inj(u(R + )) + ǫ. If n is big enough, there exists an integer k n such that B ∞ (γ -kn n i, i) belongs to I.

Proof of

Step 2. For a positive integer k we put :

r n,k := B ∞ (γ -k n i, i).
We have

r n,k = k-1 l=0 B ∞ (γ -k+l n i, γ -k+l+1 n i).
Let us prove that for n big enough, each step

B ∞ (γ -k+l n i, γ -k+l+1 n i) is smaller than Inj(u(R + )) + ǫ.
For now, we admit the following lemma :

Lemma 3.2. For all positive integers n and nonpositive integer a, there exists a point p n,a in H 2 satisfying :

(i) d(p n,a , γ n p n,a ) = d(q n , γ n q n ) (ii) B ∞ (p n,a , γ a n i) = -d(γ a n i, p n,a ) For convenience, set s n,a := d(γ a n i, p n,a ) and A k,l := B ∞ (γ -k+l n i, γ -k+l+1 n i).
Using Lemma 3.2, we compute :

A k,l = B ∞ (γ -k+l n i, γ -1 n p n,-k+l+1 ) + B ∞ (γ -1 n p n,-k+l+1 , p n,-k+l+1 ) + B ∞ (p n,-k+l+1 , γ -k+l+1 n i) ≤ d(γ -k+l n i, γ -1 n p n,-k+l+1 ) + d(γ -1 n p n,-k+l+1 , p n,-k+l+1 ) -s n,-k+l+1 = s n,-k+l+1 + d(q n , γ n q n ) -s n,-k+l+1 = d(q n , γ n q n )
As d(q n , γ n q n ) is eventually smaller than Inj(u(R + )) + ǫ, we obtain that for n big enough :

A k,l = B ∞ (γ -k+l n i, γ -k+l+1 n i) ≤ Inj(u(R + )) + ǫ.
(3) Thus, r n,k is a sum of k terms A k,l , for l = 0 . . . k -1, all smaller than Inj(u(R + )) + ǫ which is the length of the interval I of R + . We now use the fact that since γ -k n i -→ k→+∞ γ - n = ∞ for every n, we have r n,k -→ k→+∞ +∞ for every n. Thus, there exists an integer k n such that r n,kn belongs to I. This concludes Step 2.

End of the proof of Theorem 1.1. Take the elements γ n given by Step 1. Fix an ǫ > 0, chosen arbitrarily small, and an interval I of R + of length Inj(u(R + ))+ ǫ. Consider the sequence of positive integers (k n ) N given by Step 2. Since the numbers B ∞ (γ -kn n i, i) eventually all belong to I, the sequence of numbers (B ∞ (γ -kn n i, i)) N admits an accumulation point r in I. Thus, setting α n := γ kn n and applying the equivalence (1), we get that g r u ∈ h R uh R u. Now, applying the same argument to a partition (I l ) N of R + in intervals of lengths Inj(u(R + )) + ǫ, we get a sequence of times (t l ) N , where each t l is in I l , and such that g t l u belongs to h R u for every l.

Moreover, if Inj(u(R + )) = 0, as the intervals I l will be of length 0 + ǫ for an ǫ > 0 arbitrarily small, we get g R + u ⊂ h R u.

Proof of Lemma 3.2. Take an isometry γ n as in Step 1.

Observe that a point in C γn (q n ) ∩ [γ a n i, ∞) would satisfy condition (i) and also condition (ii) according to Proposition 2.1 and Proposition 2.2. We prove that this intersection is not empty.

Suppose that γ n is hyperbolic. There are two cases. The first case is when the signs of γ - n and γ + n are opposed. Since lim n→+∞ γ - n = lim n→+∞ γ + n = ∞, the graph of C γn (i) is below the one of C γn (q n ). So every geodesic starting from a point z on C γn (i) and ending at ∞ has an intersection with C γn (q n ). This is true in particular when z = γ a n i. The second case is when γ - n and γ + n have same sign. Observe then that we eventually have

d(i, (γ - n , γ + n )) > d(q n , (γ - n , γ + n ))
, because the converse would contradict the discretness of Γ using Proposition 2.1. So, according to the same proposition, the graph of C γn (q n ) is below the one of C γn (i) if n is big enough. Observe now that as a is a nonpositive integer, the point γ a n i belongs to the portion of C γn (i) between i and γ - n , and that for any point z in this portion of

C γn (i), the intersection C γn (q n ) ∩ [z, ∞) is not empty.
If γ n is parabolic, the proof is similar to the second case, when replacing

d(i, (γ - n , γ + n )) by B γ + n (i, C γn (z 0 )) and d(q n , (γ - n , γ + n )) by B γ + n (q n , C γn (z 0 )
) and using Proposition 2.2. This concludes the proof.

An example to prove theorem 1.3

Recall that the Dirichlet domain centered at i of a fuchsian group Γ with no elliptic element fixing i is defined by :

D i (Γ) := γ∈Γ γ =id H i (γ),
where

H i (γ) := {z ∈ H 2 | d(z, i) ≤ d(z, γ(i))}.
The main interest for us of looking for this domain is the following classical fact (see section I, Proposition 4.9 of [START_REF] Dal | Geodesic and Horocyclic Trajectories[END_REF]

) : if u ∈ Γ\(T 1 H 2 ) and if for some lift ũ in T 1 H 2 the point ũ+ belongs to D i (Γ) ∩ ∂H 2 , then u(R + ) is almost minimizing.
Let us now construct our example. Consider, for any rational number q ∈ [4, +∞) and any n in N, the hyperbolic isometry :

g q,n := √ q (1 -q)r n -1 rn √ q ,
where (r n ) N is the sequence of real numbers defined by :

r 1 = 2 r n = 3r n-1 , ∀n ≥ 2. Let F 1 := {g q,n , q ∈ Q ∩ [4, +∞), n ∈ N}.
We now conjugate the isometries g 4,n by the isometries T q,n := 1 t q,n 0 1 , for any rational number q ∈ (1, 4) and any n in N, with t q,n := -r n ( √ q -2). We have :

h q,n := T -1 q,n g 4,n T q,n = 4 - √ q r n [( √ q -2) 2 -3] -1 rn √ q .
Set F 2 := {h q,n , q ∈ Q ∩ (1, 4), n ∈ N}. Define, for every q ∈ Q ∩ (1, +∞), the hyperbolic isometry f q,n by :

• f q,n := h q,n if q ∈ Q ∩ (1, 4), • f q,n := g q,n if q ∈ Q ∩ [4, +∞),
and set F := F 1 ∪ F 2 = {f q,n , q ∈ Q ∩ (1, +∞)}.
For any non-elliptic isometry γ in PSL 2 (R), define ∂H i (γ) to be the perpendicular bisector of the segment [i, γi]. Also denote by c(γ) the center of the euclidean half-circle ∂H i (γ) and by e l (γ), with l = 1, 2, the extremities in ∂H 2 of ∂H i (γ). Denote also by C (η,ξ) (z) the hypercycle with extremities η and ξ in ∂H 2 and passing through z in H 2 .

The following key proposition gives us all the necessary information about the perpendicular bisectors ∂H i (γ) (see section 5 for the proof).

Proposition 4.1. For every q ∈ (1, +∞) we have the following :

(i) lim n→+∞ c(f q,n ) = -∞ and lim n→+∞ c(f -1 q,n ) = +∞. (ii) lim n→+∞ e l (f q,n ) = -∞ and lim n→+∞ e l (f -1 q,n ) = +∞ for l = 1, 2. (iii) Each perpendicular bisector ∂H i (f -1 q,n ) is below the hypercycle C (0,∞) (1+i)
. Using Proposition 4.1, our goal is to extract a sequence of elements γ m of F such that all the perpendicular bisectors ∂H i (γ m ) are disjoint and which contains an infinite number of elements f q,n for every rational number q ∈ (1, +∞).

Consider any bijection ψ : N → Q∩(1, +∞)×N and set ψ(m) = (q m , ψ 2 (m)). Observe that for every q in (1, +∞), there is an infinite number of elements q m such that q m = q. We can now write F = {f qm,ψ2(m) , m ∈ N}.

Choice of the elements γ m . Put γ 0 := f q0,0 and set 2C to be the distance between the geodesic (0, ∞) and the hypercycle C (0,∞) (1 + i). Then, for every m ≥ 1, we define γ m by induction. We ask that γ m = f qm,nm where n m is the smallest integer among the integers p satisfying :

(i) |e ǫ k (f qm,p )| > |e ǫ l (γ m-1 )|, for k = 1, 2, l = 1, 2 and ǫ = ±1. (ii) d(∂H i (f qm,p ), ∂H i (f -1 qm,p )) ≥ 2C (iii) d(∂H i (f ±1 qm,p ), ∂H i (f ±1 qs,ns )) ≥ C for all s < m.
Such a sequence (γ m ) N exists according to Proposition 4.1. We now set Γ :=< γ m , m ∈ N > and prove that Γ answers Theorem 1.3.

By a classic ping-pong argument (see [START_REF] Dal | Geodesic and Horocyclic Trajectories[END_REF]), the group Γ is discrete and free. Moreover, its Dirichlet domain centered at i is :

D i (Γ) = m∈N H i (γ m ) ∩ m∈N H i (γ -1 m ).
Fix ũ in T 1 H 2 such that ũ(0) = i and ũ+ = ∞ and consider its projection u to T 1 S := Γ\(T 1 H 2 ). Since ∞ is in D i (Γ) ∩ ∂H 2 , the half-geodesic u(R + ) is almost minimizing. Hence (i) of the theorem.

To prove condition (ii), we use Proposition 2.3. Observe that it follows directly from the definition of f q,n that for every rational number q in (1, +∞), we have

f -1 q,n ∞ -→ n→+∞ ∞. (4)
Moreover, since Im(f q,n (i)) -→ n→+∞ 1 q , we have :

B ∞ (f q,n i, i) -→ n→+∞ B ∞ ( i q , i) = | ln(q)|. (5)
Now, for every q ∈ Q ∩ (1, +∞) there is an infinite number of elements f q,n in (γ m ) N , thus in Γ. So, according to (4), (5) and Proposition 2.3, it follows that all the elements g | ln(q)| u belong to h R u. Hence, g R + u is included in h R u and we get (ii) of the theorem.

Finally, fix z in the interior of D i (Γ) and any γ = γ i1 m1 . . . γ i k m k in Γ different from the identity, written as a reduced word in the letters γ m .

If k = 1, then γ = γ i1 m1 and

d(z, γz) = d(z, γ -1 z) ≥ max(d(z, ∂H i (γ m1 )), d(z, ∂H i (γ -1 m1 ))) Since d(∂H i (γ m1 ), ∂H i (γ -1 m1 )) ≥ 2C, we obtain that d(z, γz) ≥ C. If k > 1, d(z, γz) = d(z, γ i1 m1 . . . γ i k m k z) = d(γ -i1 m1 z, γ i2 m2 . . . γ i k m k z) ≥ d(∂H i (γ ±1 m1 ), ∂H i (γ ±1 m2 )) ≥ C.
It follows that the injectivity radius on S is everywhere greater than C. Finally, since all the axes of the elements γ m in Γ intersect the half-geodesic [i, ∞), and since their translation length l γn is constant by definition of the elements f q,n , the injectivity radius Inj(u(R + )) is also finite. So C ≤ Inj(u(R + )) < +∞. Hence condition (iii) of Theorem 1.3. This completes the proof.

Proof of Proposition 4.1

Using the classic formula :

∀a, b ∈ H 2 , sinh d(a, b) 2 = |a -b| 2 Im(a)Im(b) ,
we get :

Proposition 5.1. Consider a point P = R + iI in H 2 . The equation of the perpendicular bisector of the hyperbolic segment between i and P is :

x + R I -1 2 + y 2 = I 1 + R 2 (I -1) 2
For the following calculations, we distinguish the case q ∈ [4, +∞) from the case q ∈ (1, 4).

Case 1: fix q ∈ [4, +∞). We have f q,n i = R q,n + iI q,n where

R q,n := √ q(1 -q)r n - √ q rn q + 1 r 2 n and I q,n := 1 q + 1 r 2 n .
Observe that we have r n -→ n→+∞ +∞. Thus, as n converges to +∞, the quantity R q,n is equivalent to √ q(1-q)rn q = r n (1-q) √ q , where the number 1-q √ q is different from 0, and the quantity I q,n is equivalent to 1 q . So, applying Proposition 5.1, we get the following asymptotic equivalence :

c(f q,n ) = - R q,n I q,n -1 ≍ n→+∞ -r n √ q.
So the centers c(f q,n ) converge to -∞.

Let us now study the radii of the geodesics ∂H i (f q,n ). We have

I q,n 1 + R 2 q,n (I q,n -1) 2 1 2 ≍ n→+∞ 1 √ q R q,n I q,n -1 = - 1 √ q c(f q,n )
where 1 √ q belongs to (0, 1 2 ].

Thus, according to Proposition 5.1, the extremities e l (f q,n ), for l = 1, 2, of the geodesics ∂H i (f q,n ) converge to -∞, and since 1 √ q ≤ 1 2 , all these geodesics are below the hypercycle C (0,∞) (-1 + i).

We now study the case of ∂H i (f -1 q,n ). We have :

f -1 q,n i = √ q(q -1)r n + √ q rn q + 1 r 2 n + i 1 q + 1 r 2 n .
So we observe that the real part of f -1 q,n i is the negative of the real part of f q,n i. So the geodesics ∂H i (f -1 q,n ) and ∂H i (f q,n ) are symetric with respect to the imaginary axis. In particular, they are below the hypercycle C (0,∞) (1 + i).

Case 2: fix q ∈ Q ∩ (1, 4). We have f q,n i = R q,n + iI q,n where

R q,n := r n √ q(( √ q -2) 2 -3) -1 rn (4 - √ q)
q + 1 r 2 n and I q,n = 1 q + 1 r 2 n .

Observe that since q ∈ (1, 4), the number √ q(( √ q -2) 2 -3) is different from 0. So as n goes to +∞, we have the equivalences : R q,n ≍ n→+∞ r n ( √ q -2) 2 -3 √ q and I q,n ≍ n→+∞ 1 q .

Thus, according to Proposition 5.1 c(f q,n ) = -R q,n I q,n -1 ≍ n→+∞ -r n √ q ( √ q -2) 2 -3 1q where √ q

( √ q-2) 2 -3 1-q
is a real number greater than 2. Thus, these centers converge to -∞.

Let us now study the radii of the geodesics ∂H i (f q,n ). We have

I n,q 1 + R 2 n,q (I n,q -1) 2 1 2 ≍ n→+∞ 1 √ q R n,q I n,q -1 = - 1 √ q c(f q,n )
where 1 √ q belongs to ( 1 2 , 1). Thus, again, the extremities e l (f q,n ), for l = 1, 2, of the geodesics ∂H i (f q,n ) converge to -∞ as n goes to +∞.

We now study the case of f -1 q,n for q ∈ (1, 4). We have f -1 q,n i = R q,n + iI q,n where R q,n = -r n (4 -√ q)[( √ q -2) 2 -3] + √ q rn (4 -√ q) 2 + 1 r 2 n and I q,n = 1 (4 -√ q) 2 + 1 r 2 n .

Since q ∈ (1, 4), the number (4 -√ q)[( √ q -2) 2 -3] is different from 0. Thus,

R q,n ≍ n→+∞ r n ( √ q -2) 2 -3 √ q -4 . Since I q,n ≍ n→+∞ 1 (4 - √ q) 2 ,
we have :

c(f -1 q,n ) = - R q,n I q,n -1 ≍ n→+∞ -r n (( √ q -2) 2 -3)( √ q -4) 1 -(4 - √ q) 2
Since the number

(( √ q-2) 2 -3)( √ q-4) 1-(4- √ q) 2
is negative for q ∈ (1, 4), the centers c(f -1 q,n ) of the geodesics ∂H i (f -1 q,n ) converge to +∞ as n goes to +∞. We now study the radii :

I n,q 1 + R 2
n,q (I n,q -1) 2 1 2 ≍ n→+∞ 1 4 -√ q R q,n I q,n -1 = -1 4 -√ q c(f -1 q,n )

where the number 1 4-√ q belongs to ( 1 3 , 1 2 ). Thus, the extremities e l (f -1 q,n ), for l = 1, 2, converge to +∞. Moreover, since 1 4-√ q < 1 2 , the geodesics ∂H i (f -1 q,n ) are below the hypercycle C (0,∞) (1 + i) as claimed.