Experimental and simulation-based investigation of the influence of freezing and annealing on the microstructure in lyophilisates

EXPERIMENTAL AND SIMULATION-BASED INVESTIGATION OF THE INFLUENCE OF FREEZING AND ANNEALING ON THE MICROSTRUCTURE IN LYOPHILISATES Doctoral Thesis Submitted in partial fulfilment of the requirements for the degree of doctor of natural sciences of the Faculty of Mathematics and Natural Sciences at Kiel University, Germany by Tigran Levonovic Kharatyan Kiel, 2024 Printed with the permission of the Faculty of Mathematics and Natural Sciences of Kiel University. This work was carried out from May 2019 until April 2023 at the Department of Pharmaceutical Development of Daiichi Sankyo Europe GmbH, Pfaffenhofen an der Ilm, Germany. The thesis was prepared under supervision of Prof. Dr. Regina Scherließ. In this thesis, DeepL Write has been used to identify and correct language mistakes (grammar and spelling). 1st Referee: Prof. Dr. Regina Scherließ 2nd Referee: Dr. Nora A. Urbanetz Submission of the PhD application: 21.05.2024 Date of examination: 06.09.2024 Dean: Prof. Dr. Frank Kempken – Dedicated to my family – Research activities Publications The dissertation includes text and figures from the publications of the authors listed below. I. Quantitative Analysis of Glassy State Relaxation and Ostwald Ripening during Annealing Using Freeze-Drying Microscopy by Tigran Kharatyan, Srikanth R. Gopireddy, Toru Ogawa, Tatsuhiro Kodama, Norihiro Nishimoto, Sayaka Osada, Regina Scherließ and Nora A. Urbanetz in Pharmaceutics (Volume 14, Issue 6, June 2022, 1176) II. Impact of post-freeze annealing on shrinkage of sucrose and trehalose lyophilisates by Tigran Kharatyan, Shunya Igawa, Srikanth R. Gopireddy, Toru Ogawa, Tatsuhiro Kodama, Regina Scherließ and Nora A. Urbanetz in International Journal of Pharmaceutics (Volume 641, June 2023, 123051) Table of contents Table of contents Introduction ..................................................................................................................... 1 Objectives......................................................................................................................... 2 1 General background ....................................................................................................... 3 1.1 Freeze-drying in the pharmaceutical field ........................................................................ 3 1.2 Freeze-dryer unit .............................................................................................................. 4 1.3 Pharmaceutical formulations in freeze-drying ................................................................. 5 1.4 Process steps in a freeze-drying cycle .............................................................................. 7 1.4.1 Freezing ...................................................................................................................... 9 1.4.1.1 Heterogeneous and homogeneous nucleation ................................................. 10 1.4.1.2 Freeze-concentration of the matrix phase........................................................ 12 1.4.2 Annealing during the freezing phase ....................................................................... 15 1.4.3 Primary drying .......................................................................................................... 17 1.4.4 Secondary drying ..................................................................................................... 20 1.5 Microstructure formation in an aqueous solution during freezing................................ 21 1.6 Coarsening of the ice crystal phase during annealing .................................................... 23 1.7 Morphological assessment of the microstructure in a lyophilisate ............................... 25 1.8 Numerical simulations and phase-field methods in the pharmaceutical field .............. 27 2 Materials and methods ................................................................................................ 31 2.1 Materials ......................................................................................................................... 31 2.1.1 D-(+)-sucrose ............................................................................................................ 31 2.1.2 D-(+)-trehalose ......................................................................................................... 32 2.2 Methods .......................................................................................................................... 32 2.2.1 Lyophilisation ........................................................................................................... 32 2.2.2 Freeze-drying microscopy ........................................................................................ 33 2.2.3 Polarised light microscopy ....................................................................................... 35 2.2.4 Differential scanning calorimetry ............................................................................ 36 2.2.5 Phase-field computation and visualisation .............................................................. 36 2.2.6 Image analysis .......................................................................................................... 37 3 Results and discussion .................................................................................................. 38 3.1 Microstructure formation during freezing ..................................................................... 38 Table of contents 3.1.1 Recalescence and crystallisation in samples during freezing .................................. 38 3.1.2 Characterisation of the thermal profile during freezing.......................................... 41 3.1.3 Height-dependent recalescence and crystallisation in lyophilisation vials ............. 48 3.1.4 Assessment of the microstructure in a frozen solution ........................................... 52 3.1.5 Impact of thermal history during freezing on the microstructure .......................... 62 3.1.6 General discussion – importance of recalescence and crystallisation .................... 66 3.2 Microstructure evolution during annealing ................................................................... 67 3.2.1 Methodological framework to investigate the microstructure ............................... 68 3.2.2 Implementation and calibration of the phase-field model ..................................... 72 3.2.2.1 Formulism of the multiphase-field model ........................................................ 72 3.2.2.2 Experimental determination of relevant calibration parameters..................... 74 3.2.2.2.1 Mass fractions during annealing ............................................................... 75 3.2.2.2.2 Volume fractions during annealing/maximum freeze-concentration ...... 78 3.2.2.2.3 Recrystallisation rates during annealing ................................................... 80 3.2.3 Initial conditions for two-dimensional and three-dimensional simulations ........... 86 3.2.4 Calibration of input parameters for the phase-field simulations ............................ 95 3.2.4.1 Adjustment of the area and volume fractions of the particulate phase .......... 95 3.2.4.2 Coarsening parameters for the two-dimensional phase-field model ............. 101 3.2.5 Advancement of the phase-field simulation to three-dimensional ...................... 107 3.2.6 Final adjustment of particle volume fractions ....................................................... 110 3.2.7 Morphological properties in the three-dimensional phase-field model ............... 111 3.2.7.1 Description of morphological properties in a porous microstructure ............ 112 3.2.7.2 Ice crystal size distribution .............................................................................. 116 3.2.7.3 Conjunction area of ice crystals ...................................................................... 120 3.2.7.4 Surface area of the matrix phase .................................................................... 126 3.2.8 General discussion – pore size determination in lyophilisates ............................. 129 4 Future applications ..................................................................................................... 134 4.1 Reorganisation kinetics during recalescence and crystallisation ................................. 134 4.2 CFD simulations through the porous microstructure ................................................... 135 5 Appendices ................................................................................................................. 137 Appendix A – Validation and stability analysis for the phase-field simulation .................. 137 i. Gibbs-Thomson (capillarity) effect on composition .................................................... 137 ii. Energy barrier coefficient on interfacial free energy .................................................. 138 General background 5 The solution intended for drying is placed in containers such as glass vials, which are available in various sizes and materials. For instance, amber glass vials are used for products sensitive to UV radiation [17]. Before drying, the product must be thoroughly frozen at temperatures typically ranging from -50 °C to -80 °C. This ensures maximum phase separation within the sample due to freezing, a process that will be detailed later. Consequently, a cooling system connected to the shelves where the vials are placed is essential. Some freeze-dryers lack an integrated cooling system, requiring products to be pre-frozen outside the dryer. To facilitate sublimation, the pressure within the drying chamber must be reduced to at least below the triple point of the solvent, making a vacuum system indispensable. The drying chamber is connected to a condenser chamber and a vacuum pump through valves, with the vacuum pump regulating the pressure by removing non-condensable gases [16]. The condenser, kept at temperatures lower than the shelves, enables the transition of water vapor into ice via deposition/resublimation. The driving force behind the freeze-drying process is the temperature difference between the product and the condenser. This difference creates varied saturation vapor pressures of ice at the sublimation front in the product and at the ice on the condenser, thereby establishing a pressure gradient that facilitates the transfer of water molecules from the product to the condenser [18]. Additionally, these condensers are designed with high surface areas to enhance their capacity for accumulating ice [2]. The vials are heated via the shelves in order to expedite the sublimation process by further increasing the temperature difference between the product and the condenser [19]. 1.3 Pharmaceutical formulations in freeze-drying Freeze-drying is typically utilised for Active Pharmaceutical Ingredients (APIs) where gentle drying processes and stabilisation are crucial [3]. This is particularly the case for biopharmaceutical APIs such as monoclonal antibodies, therapeutic proteins and vaccines, which are prone to loss of functionality, expensive to produce, and require meticulous downstream processing [5,20]. Furthermore, temperature-sensitive small molecule drugs such as antibiotics or antivirals are also subject to freeze-drying in order to preserve their functionality. Some examples of freeze-dried pharmaceuticals are listed in Table 1.1. Besides the API, freeze-dried pharmaceutical products contain inactive excipients, defined by the FDA as “any component of a drug product other than an active ingredient” [21]. General background 6 Approximately 67 % of lyophilised small molecule pharmaceuticals and nearly all biotechnology-derived products contain excipients [22,23]. Lyophilisation excipients are classified according to their functions as stabilisers, buffers, surfactants, or bulking agents [4]. Certain substances may fulfil multiple roles in a formulation. For example, disaccharides such as sucrose may be used as stabilisers and bulking agents at the same time [24]. Table 1.1: Examples of freeze-dried pharmaceuticals with FDA approval [25]. Drug name (FDA approval year) API Category/Indication Manufacturer Enbrel® (1998) Etanercept Rheumatoid arthritis, Psoriasis Amgen Inc. Enhertu® (2019) Trastuzumab deruxtecan Breast/gastric cancer Daiichi Sankyo Company, Limited Avonex® (1996) Interferon beta-1a Multiple sclerosis Biogen Orencia® (2005) Abatacept Rheumatoid arthritis Bristol-Myers Squibb Botox® (1991) Daxibotulinumtoxin A Various Allergan Dantrium® (1979) Dantrolene sodium Muscle relaxant The Procter & Gamble Company Cosmegen® (2009) Dactinomycin Antibiotic oncologic Merck KGaA Stabilisers serve to protect APIs from aggregation or structural changes during freezing or drying, acting as cryoprotectants or lyoprotectants [4,26–28]. Examples of stabilisers include disaccharides such as sucrose and trehalose, along with synthetic polyethers such as polyethylene glycol and polyvinylpyrrolidone. The mechanisms through which stabilisers operate include vitrification and water replacement [29]. Vitrification involves the immobilisation of the API in an amorphous matrix phase, which prevents interactions between API molecules. Water replacement describes the substitution of hydrogen bonds between the API and water with bonds between the stabiliser and the protein. As the stabiliser is not removed during drying, the API is less likely to be forced out of its native form upon desorption of unfrozen water. Buffers can be incorporated into a formulation to regulate pH and maintain the solubility and stability of the API. For instance, proteins typically have the lowest solubility at their isoelectric point due to the neutral net charge on their surface, which leads to interactions between General background 7 proteins rather than between proteins and dipolar water molecules [30]. However, the selection of a suitable buffer is crucial, as phosphate buffers, for example, can lead to significant pH shifts during freezing, resulting in protein denaturation and precipitation [31,32]. Buffers commonly used in lyophilisation, such as citrates and histidine buffers, are less prone to crystallisation and exhibit a smaller temperature-dependent pH range [4]. Surfactants, such as polysorbate 80 and lecithin, are amphiphilic and effective at lowering the surface tension of water. In pharmaceutical formulations, surfactants bind to the hydrophobic regions of proteins, decreasing aggregation. Furthermore, surfactants raise the energy required for proteins to unfold, thereby stabilising their native structure [33,34]. Bulking agents are added to formulations when the API is present in small amounts per dose due to its high potency, therefore increasing the total mass of the pharmaceutical product [35]. Frequently utilised bulking agents for biopharmaceuticals include disaccharides like sucrose and trehalose, or sugar alcohols such as xylitol and mannitol, which form either amorphous or crystalline phases upon freezing [2,36,37]. As bulking agents may constitute a significant portion of a lyophilisate, their behaviour during freezing and drying is important for the overall process. 1.4 Process steps in a freeze-drying cycle A freeze-drying cycle is divided into three main process steps: freezing, primary drying, and secondary drying. During each stage, the pharmaceutical solution or lyophilisate undergoes physicochemical changes on both the microand macroscale, as illustrated in Figure 1.3. The subprocesses of a lyophilisation cycle are regulated by the shelf temperature and chamber pressure. The initial freezing step involves lowering the temperature of the samples until the solution solidifies completely, resulting in the formation of a microstructure composed of dispersed individual ice crystals within a continuous matrix phase. Solutions containing glassforming solutes freeze over a wide temperature range. For example, a solution might start as a slush around -10 °C and solidify around -30 °C to -40 °C, depending on the solute [38]. Therefore, to thoroughly freeze such a solution, temperatures well below the triple point of the solvent are necessary. An optional annealing subprocess can be included in the freezing step, where the temperature is slightly raised but kept below the melting point of the solution to prevent thawing. General background 8 Figure 1.3: Illustration of a freeze-drying process including an annealing step during freezing. This phase is shown in greater detail to highlight the annealing temperature (Ta), annealing time (ta), and the rates of freezing and cooling. In primary drying, most of the water, in the form of ice crystals, is removed from the vials by reducing the chamber pressure to at least below the triple point of the solvent. Since sublimation is an endothermic phase transition, the temperature inside the vials decreases during this stage. To counteract this cooling and accelerate sublimation, the samples are heated via the shelves. During this phase, only the ice in direct contact with the gas phase undergoes the phase transition, creating a sublimation front that progresses from the surface of the frozen solution towards the bottom of the vial. To complete the drying process, it is necessary to remove the residual water from the matrix phase. This is typically achieved by further increasing the shelf temperature, which causes water molecules to diffuse and desorb from the matrix phase, reducing the water content in the lyophilisate to the target range of usually around 1 % (w/w) to 3 % (w/w) for freeze-dried biological products [2]. Once the desired water content is reached, the vials are sealed, often with an inert gas, to prevent moisture from the air from wetting the lyophilisate. General background 9 1.4.1 Freezing In the freezing step, the shelf temperature is lowered with the goal of ensuring that most of the solvent, such as water, crystallises. This leads to partial phase separation, commonly referred to as freeze-concentration, which is crucial for the subsequent removal of the solvent through sublimation in later stages of the process. This phase includes several subprocesses, which are illustrated in Figure 1.4. Figure 1.4: Temperature profile of sample and shelf during the freezing step: precooling (1), freezing/crystallisation (2), freezing plateau (3), annealing ramp (4), annealing (5), refreezing (6), refreezing plateau (7). After an optional precooling step aimed at equalizing the conditions across vials in the freezedryer to approximately 5 °C, the shelf temperature is further decreased to typically between -40 °C and -50 °C to freeze the samples. During this process, a notable increase in the sample's temperature occurs after reaching 10 °C to 15 °C below its equilibrium freezing/melting point, i.e., the temperature at which a frozen solution with the same components would completely melt. This temperature increase signals an exothermic phase transition, specifically the crystallisation of water to ice. Just before the temperature begins to rise, the solution remains liquid despite being below its equilibrium freezing temperature, thereby becoming supercooled. A supercooled solution is not in equilibrium and is only meta-stable [39]. The greater the degree of supercooling, the more the solution becomes oversaturated with respect to the solvent, which in turn increases the probability of nucleation [13]. Once General background 10 nucleation begins, exothermic crystallisation rapidly occurs throughout the entire sample. This rise in temperature during solidification is known as recalescence, occurring when the heat released during the transition exceeds the rate at which heat can dissipate from the material [40]. This phenomenon, depicted in Figure 1.5, is essential for understanding the impact of freezing conditions on the lyophilisation process. Figure 1.5: Schematic temperature profile of a solution during freezing: supercooling (1), nucleation (2), recalescence (3), freeze-concentration (isothermal crystallisation) (4), end of solidification (5). The equilibrium freezing temperature (Tf) is shown as a grey dashed line. Modified from [41]. Once the sample temperature nearly reaches its equilibrium freezing temperature, it remains steady for period of time. Meanwhile, water continuously crystallises, releasing heat that must be removed from the sample by the cooling shelves. A detailed description of the phenomena contributing to the typical temperature profile of the sample during freezing will be presented in the following sections. Additionally, an in-depth examination of recalescence and the subsequent crystallisation will be explored in the Results section of this work. This focus is motivated by the often underestimated impact of the temperature increase during freezing on the microstructure of the lyophilisate. 1.4.1.1 Heterogeneous and homogeneous nucleation The onset of crystallisation requires the formation of initial nucleation points, known as nuclei, where water molecules start to precipitate and form larger ice crystals [42]. Nucleation does General background 11 not necessarily occur at the equilibrium freezing temperature but can happen at much lower temperatures, influenced by process conditions and material properties. For instance, pure water, with a melting point at 0 °C, can remain unfrozen at temperatures as low as -42 °C [43]. However, the presence of impurities might initiate nucleation at temperatures closer to the equilibrium freezing point, such as -5 °C [44,45]. These impurities could include particulate contaminants or large molecules like proteins, acting as interfaces to facilitate nucleation [46]. This phenomenon is explained by classical nucleation theory, which differentiates between homogeneous and heterogeneous nucleation, and is illustrated in Figure 1.6. Figure 1.6: Size-dependent free energy of homogeneous (A) and heterogeneous (B) nucleation. Modified from [47]. In homogeneous nucleation, small nuclei form as a solid phase within the oversaturated solution. These nuclei are clusters of water molecules arranged and bonded similarly to ice. Their formation, driven by Brownian motion within the bulk volume, is stochastic in nature, yet the probability of formation increases as the temperature decreases [48]. The energy required to form and sustain a nucleus depends on its size, with the Gibbs free energy (ΔG) comprising surface free energy and bulk free energy components. Surface free energy increases with nucleus size, while bulk free energy has a stabilising effect, reducing the system's free energy as the nucleus grows [49]. Therefore, the stability of a nucleus hinges on its surface area to volume ratio, or simply its radius in the case of a spherical nucleus. The free energy peak indicates the critical particle size (rc) that must be surpassed for larger ice crystals to grow from the nucleus. Failing to exceed this size, the system's free energy can only decrease by reducing the surface area of the nuclei, potentially leading to their complete break up. In heterogeneous nucleation, the new solid phase forms on an existing interface, General background 12 such as the surface of an impurity, which reduces the required energy for nucleus stabilisation due to a smaller surface area compared to homogeneous nucleation for the same radius, thereby making heterogeneous nucleation more prevalent [50]. In freeze-drying, various techniques are utilised to influence the nucleation temperature during the freezing step. These include the ice fog technique, quench freezing, electro freezing, ultrasound-controlled ice nucleation, vacuum-induced surface freezing, high-pressure shift and depressurisation, and the addition of nucleation agents [44]. 1.4.1.2 Freeze-concentration of the matrix phase Once stable nuclei with radii above the critical size are formed, water molecules from the supercooled liquid phase can attach to their surfaces and start to form ice crystals [51]. In the conditions of a freeze-dryer, ice typically forms in the Ih phase, also known as normal hexagonal crystalline ice. Within a solution, the dense structure of hexagonal ice prevents most solutes from being incorporated into the crystalline lattice, leading to the formation of a nearly pure ice phase and a solute-enriched matrix phase [16]. This reversible process, known as freeze-concentration, is governed solely by temperature. The lower the temperature, the more water precipitates from the solution into the existing ice crystals, as illustrated in Figure 1.7. However, this process can only continue until the matrix phase reaches its maximum freeze-concentrated state, which varies depending on the solute in the aqueous solution [52]. Figure 1.7: Freeze-drying microscopy of a partially frozen 10 % (w/w) sucrose solution at -4 °C (left) and at -10 °C (right). The freezing behaviour of a matrix phase can be classified as crystallising or vitrifying and the choice of excipients in the formulation will determine which of these behaviours is present or General background 13 whether a combination of both is observed. Crystalline systems are typically characterised by the formation of a eutectic mixture upon freezing, as illustrated in the phase diagram for a two-component system shown in Figure 1.8. Figure 1.8: Phase diagram of a binary system with crystallising components A and B. Modified from [16]. A eutectic mixture has a lower freezing point than any other mixing ratio of its components. This means that the mixture abruptly freezes at the eutectic temperature (Teu), and the components, A and B, form a fine structure of individual crystalline phases, illustrated in Figure 1.8 as α and β. Non-eutectic mixtures, either hypoeutectic or hypereutectic, begin to form a precipitate of the excess component as the temperature decreases. A hypoeutectic mixture remains entirely liquid until it reaches the equilibrium freezing temperature, also known as the liquidus line. At this point, component A begins to precipitate (1→2). This process continues along the liquidus line until the eutectic point (5), where a portion of component A forms a eutectic mixture with all of component B. Upon further cooling, this eutectic mixture crystallises. A hypereutectic mixture follows a similar path (3→4→5). Upon solidification, such eutectic systems may form microscopic arrangements, where one phase is embedded as lamellar or globular structures within the other phase [53,54]. Many pharmaceutical formulations contain major components that do not crystallise upon freezing but instead remain in a rubbery phase and undergo glass transition upon sufficient cooling. During freezing of a solution with a glass-forming component, precipitation of water General background 14 and ice formation follow the liquidus line once nucleation occurs, similar to a crystallising system. However, instead of forming a eutectic mixture, i.e., a mixture with a single melting point, more water crystallises over a broad temperature range. This process continues until the amorphous phase cannot be further concentrated by temperature reduction. At this point, the system is referred to as maximally freeze-concentrated and is characterised by the onset temperature for ice melting Tm’ and the solute concentration in the matrix phase cg’ [38]. The phase diagram in Figure 1.9 illustrates the freezing behaviour of such solutions. Figure 1.9: Phase diagram of an aqueous system with a glass-forming solute. The glass transition temperature of the solute-enriched phase (Tg) changes during freeze-concentration (blue) and desorption of unfrozen water (red). Modified from [55,56]. In case of glass-forming disaccharides, the formation of hydrogen bonds between solute and water as well as between solute molecules exhibits comparable energy and configuration [57]. A decrease in temperature shifts the chemical potential (µ) of these hydrogen bonds. Once the matrix phase becomes maximally freeze-concentrated, the µ of hydrogen bonds between the solute and water falls below the µ of water-water interactions, thus inhibiting further crystallisation of water [38]. Roos (1991) showed that approximately 20 % (w/w) of residual water remains in a maximally freeze-concentrated matrix of a frozen sucrose solution. Further investigation by Seifert et al. (2020) determined that the residual water content in sucrose ranges between 22.6 % and 24.6 % (w/w). General background 21 freeze-concentrated matrix before drying [67,94]. The conditions and duration of secondary drying allow the desired residual moisture in the lyophilisate to be achieved. This is important because APIs are not necessarily best stabilised when the residual moisture is as minimal as possible [95–97]. 1.5 Microstructure formation in an aqueous solution during freezing The literature on the properties of the porous microstructure in lyophilisates generally agrees that pore size and number are influenced by two key physical phenomena: nucleation and Ostwald ripening. The nucleation temperature is crucial as it affects the number of stable nuclei that form and subsequently grow into ice crystals. A lower nucleation temperature results in a greater number of stable nuclei [13,16,86,98]. This is further supported by molecular simulations, which indicate that for homogeneous ice nucleation at 15 °C below the equilibrium freezing temperature, the critical size of a nucleus is approximately 8000 molecules, corresponding to a critical radius of 4 nm. This critical size decreases to 1.7 nm with a greater degree of supercooling of 35 °C [99]. The volumetric nucleation rate B0 describes the number of stable nuclei that form in a given volume per unit time and is defined by Colucci et al. (2020) as: 𝐵  = 𝑘   𝑇  − 𝑇    , (Eq. 1) where kb and b are kinetic parameters, Tf is the equilibrium freezing temperature, and Tn is the nucleation temperature. The kinetic parameters in this model can be adjusted over a range of values until an acceptable fit is achieved between the model predictions and experimental data. The experimental data for this purpose can be determined by examining images of the dried lyophilisate with a scanning electron microscope (SEM) and quantifying the number of pores. According to Equation 1, the critical factor determining the number of stable nuclei is the difference between the equilibrium freezing temperature of the mixture and the actual nucleation temperature. However, this does not imply that nuclei form uniformly across the entire sample volume at the moment of nucleation. Instead, nucleation typically begins with the formation of at least one primary nucleus, likely through heterogeneous nucleation. From this initial site, the ice crystal front advances, with new nuclei forming via secondary General background 22 nucleation, often simply referred to as crystallisation, throughout the entire supercooled volume of the sample [86]. The rate of crystallisation, or more specifically, the secondary nucleation front velocity or the propagation rate of the ice front, is influenced by the degree of supercooling with greater supercooling resulting in a faster propagation rate [2,100]. Empirical equations have been formulated to depict the relationship between the propagation rate of the ice front and the nucleation temperature. For instance, at an average nucleation temperature of -10 °C ± 3 °C for water for injection, the propagation rate of the ice front is approximately 5.2 cm/s. However, this rate significantly increases to about 23.7 cm/s at a lower nucleation temperature of -20 °C [16]. When a solute is present, the propagation rate of the ice front slows down to the order of mm/s and varies depending on the solute concentration in the solution [101,102]. Once primary and secondary nuclei have formed, water molecules precipitate from the supercooled solution into these stable nuclei, forming ice crystals. Although this process does not change the number of crystals, it does increase their size. These ice crystals have crystallographic orientations that prevent them from coalescing when there is a mismatch between the orientations of neighbouring crystals [103]. The volume of water that precipitates is solely dependent on the ambient temperature, and this process continues until the maximum freeze-concentrated composition of the solute-enriched matrix phase is achieved [59]. Assuming that no additional nucleation occurs once the stable primary and secondary nuclei are present, the size of the resulting ice crystals primarily depends on the initial number of these nuclei, as illustrated in Figure 1.11. This rationale also underpins the methodology used to determine the parameters for Equation 1, where each individual pore is attributed to a corresponding individual ice crystal and, by extension, to a preceding individual stable nucleus. In their study, Colucci et al. (2020) measured pore sizes in a freeze-dried 5 % (w/w) sucrose solution using scanning electron microscopy (SEM), finding sizes ranging from approximately 25 µm to 50 µm based on axial position within the lyophilisate. Thomik et al. (2022) reported pore sizes of 23.73 µm ± 11.13 µm in a similar 5 % (w/w) sucrose solution, measured using micro-computed tomography (µ-CT). Additionally, Fang et al. (2020) employed low-pressure mercury intrusion porosimetry and Brunauer-Emmet-Teller (BET) analysis to determine General background 23 average pore sizes ranging from 20 µm to 60 µm in their study of a similar solution. Although the freezing protocols and equipment varied across these studies, the results consistently showed a similar magnitude with the minimum pore size around 20 µm. Figure 1.11: Schematic depiction of the influence of supercooling on the size of ice crystals in a frozen solution. These findings assume that size and number of pores in non-annealed samples are solely influenced by nucleation kinetics when no additional annealing step is included. However, this explanation simplifies the formation of microstructures in aqueous solutions during freezedrying. The phase transition during freezing generates latent heat of crystallisation, which raises the sample temperature and results in a complex thermal history throughout the freezing stage of a lyophilisation cycle. In summary, it is assumed that the number of stable nuclei formed during freezing dictates the number of individual ice crystals in the frozen solution, and consequently, the number of pores in the dried lyophilisate. This work, however, seeks to clarify the impact of the latent heat of crystallisation during recalescence and the subsequent crystallisation phase on the microstructural formation in solutions intended for freeze-drying. 1.6 Coarsening of the ice crystal phase during annealing An annealing step can be incorporated into the process post-freezing to coarsen the disperse ice phase, as depicted in Figure 1.12. Nakagawa et al. (2018) demonstrated that a 20 % (w/w) General background 24 sucrose solution, initially with an average pore size of approximately 60 µm after 1 hour of annealing at -5 °C, expanded to about 150 µm after 6 hours of annealing. Similarly, in the study by Thomik et al. (2022), a 5 % (w/w) sucrose solution annealed for 11 hours at -5 °C showed an increase in pore size from 23.73 µm ± 11.13 µm to 33.04 µm ± 26.95 µm. Figure 1.12: Schematic depiction of the influence of annealing on the size of ice crystals in a frozen solution. As described in Chapter 1.4.2, the phenomenon that occurs during annealing is called Ostwald ripening, and it is driven by the reduction of the system’s free energy. This effect is mathematically expressed by the classical Lifshitz-Slyozov-Wagner (LSW) theory of Ostwald ripening, which describes the coarsening behaviour of particles (here: ice crystals) or droplets in a matrix phase or solution. It is often used in the context of material science, especially for metal alloys, when characterising the growth of the average particle size over time. In the work of Niethammer (2008), a detailed description of assumptions and limitations regarding the LSW theory are given. The LSW theory assumes that the particles are spherical and exhibit initially a uniform size distribution. Additionally, the distance between particles is considered large compared to the individual particle size, so that diffusion of molecules from the particulate phase is the dominant mechanism for particle growth. Furthermore, the concentration in the surrounding medium or matrix phase is assumed to be constant, i.e., no concentration gradients are present in the matrix phase. The mathematical expression for the temporal change of particle radius r due to dissolution and redeposition is given by [105]: 𝑑𝑟 𝑑𝑡 = 𝐷 𝑟 󰇡 ∆ − 𝜀 𝑟 󰇢 , (Eq. 2) where D is a system-dependent diffusion constant, Δ is the difference between concentration at the particle surface and the equilibrium concentration in the liquid phase, ε is a term General background 25 containing surface tension, atomic volume of the solute, solubility, the ideal gas constant, and temperature. Under the assumption that diffusion of molecules through the liquid phase is the rate limiting factor, the following equation describing the change of the mean particle size can be derived [106]: 𝑟  ( 𝑡 )  = 𝑟   + 𝑘𝑡 , (Eq. 3) where r is the mean particle radius, r0 is the initial particle size, and k is the isothermal recrystallisation rate, given by [105]: 𝑘 = 8 𝜎 𝑐  𝑣  𝐷 9 𝑅𝑇 , (Eq. 4) where σ is the surface tension, c∞ is the solubility of the particle material, ν is the molar volume of the particle material, D is the diffusion coefficient of the particle material, R is the ideal gas constant, and T is the absolute temperature. According to Equation 3, the mean particle size increases with time, and the rate of increase is proportional to the cube of the particle radius. In an ideal system that fulfils all the requirements of LSW theory, Equation 4 can be used to calculate the recrystallisation rate, i.e., the rate at which the particulate phase in the system coarsens, from first principles. In such a case, the rate of Ostwald ripening can be deduced from the physical and chemical properties of the system. 1.7 Morphological assessment of the microstructure in a lyophilisate Various techniques are employed in the literature to assess the impact of annealing in a frozen disaccharide solution or the dried lyophilisate. Lyophilisates shrink due to desorption of the unfrozen water content in the maximum freeze-concentrated matrix phase during primary and secondary drying [91,67]. Although shrinkage primarily reflects a reduction in the bulk volume of the lyophilisate, it also implies changes in the internal microstructure. In the simplest scenario, a proportional reduction occurs in the individual structures, i.e., if the bulk volume decreases by approximately 20 %, the individual pores are similarly reduced in their volume. In other scenarios, internal structures may collapse to compensate for the overall reduction in bulk volume. Therefore, it is crucial to specify whether the discussion pertains to ice crystals in the frozen solution or pores in the dried lyophilisate. General background 26 Freeze-drying microscopy is commonly used to determine the collapse temperature of a solution during the primary drying step and was used in this work to investigate the impact of freezing and annealing on the initial formation and the subsequent evolution of the microstructure in frozen disaccharide solutions. Nakagawa et al. (2018) used a different approach to investigate the microstructure, where cross-sections of frozen rhodamine dyestained solutions were cut using a microtome, followed by measurement of ice crystal size with a light microscope. Alternatively, scanning electron microscopy (SEM) can be utilised to analyse the microstructure after drying [107,108]. Advancements in this field include scanning electron cryomicroscopy (cryo-SEM), which enables direct observation of frozen and undried samples [109,110]. Additionally, X-ray microtomography was employed by Thomik et al. (2022) to assess the pore size distribution in freeze-dried samples. However, each of these methods presents certain drawbacks when characterising and quantifying the microstructure that forms as a result of freezing and evolves during annealing. Freeze-drying microscopy, for instance, confines the sample between two glass plates rather than in a larger bulk volume, which can affect the coarsening of the ice crystal phase. Techniques requiring the slicing of samples, such as microtome use, yield results that depend on the fracture or slice plane in the frozen solution or lyophilisate. Since ice crystals or pores are exposed at random and varying positions, differentiation becomes challenging if only the narrowest or widest parts of an ice crystal or pore are visible, making all measured parameters imprecise. Similarly, SEM and cryo-SEM analyses produce two-dimensional data from a threedimensional structure due to the opacity of the sample surface, limiting observations to an area rather than a volume. Techniques that require dried samples also pose problems, as the removal of water may alter the structure within the matrix phase. Additionally, conducting an entire lyophilisation cycle just to test a single annealing condition can be excessively timeconsuming. Moreover, the resolution of methods like X-ray microtomography may not be sufficient to capture all fine structures within the lyophilisate. This limitation was illustrated by Thomik et al. (2022), who reported varying porosities for annealed versus non-annealed samples. This finding is contrary to expectations since annealing should increase the average pore size, yet the total volume of the matrix phase, and consequently the porosity, should remain largely unchanged. General background 27 Due to the limitations mentioned above, current methods for determining the microstructure in a frozen solution or in a dried lyophilisate are insufficient to conclusively assess the impact of freezing and annealing conditions on the microstructure. Consequently, a different approach is necessary. This approach should not be limited to experimental methods alone but should also incorporate computational techniques. These computational methods can simulate phenomena that are not directly observable or measurable, providing a more comprehensive understanding of the processes involved. 1.8 Numerical simulations and phase-field methods in the pharmaceutical field Numerical simulations are computational techniques that are used to model phenomena or entire systems as they evolve over a period of time. By solving algorithms and equations on a computer, these simulations make predictions about the system of interest [111,112]. These in silico methods can be used complementary to experimental findings, i.e., experimental data can be used to calibrate simulations, but processes can also be simulated when the acquisition of experimental data is not feasible. Overall, simulation techniques can provide valuable insights into complex systems and phenomena, but their meaningful use depends on appropriate consideration of limitations and assumptions [113,114]. In general, simulation models can be classified into two types: first-principle and phenomenological models. First-principle models, also known as ab initio (“from the beginning") models, are built directly on fundamental physical or chemical laws and principles. This means that the underlying equations represent the most basic, irreducible elements of these laws. The advantage of first-principle approaches is their utility in exploring areas with limited data availability, providing insights based on theoretical foundations. On the other hand, phenomenological models are empirical, developed based on observations and experimental data. These models mathematically describe phenomena, such as the rate of change of a parameter, and are calibrated to replicate observed behaviour as accurately as possible. While typically simpler than first-principle models, phenomenological models rely more heavily on assumptions, which may affect the reliability of their predictions. Various simulation techniques exist, each with its advantages and disadvantages depending on the application area and the system under investigation. These simulations can model physical, chemical, biological, and other types of systems, and are extensively used across General background 28 scientific research, engineering design, and numerous other fields. For instance, computational fluid dynamics (CFD) is widely employed to simulate the flow behaviour of fluids in three-dimensional spaces. CFD simulations are grounded in the fundamental laws of physics, as they derive from the conservation laws of mass, momentum, and energy, and do not depend on empirical relationships or assumptions. Applications of CFD in the field of lyophilisation include optimising the design of freeze-dryers [115] or condensers [116], as well as investigating the heat transfer within freeze-dryers [117,118], or the mass flux in vials [119]. Furthermore, the scale-up of processes is also supported by CFD analysis [120]. Simulation methods other than three-dimensional modeling are often employed to describe parts of the freeze-drying process. For instance, Colucci et al. (2020) developed a onedimensional population balance model that predicts the pore size in a lyophilisate based on its axial position within the vial. Similarly, Chun et al. (2020) proposed a mathematical model that correlates primary drying times with the size of ice crystals in frozen solutions. Additionally, Arsiccio et al. (2017) utilised a mechanistic model to predict the distribution of ice crystal sizes after freezing. While detailed information on three-dimensional simulations of ice crystals in pharmaceutical solutions during annealing is scarce, the formation of similar microstructures, which consist of particulate and matrix phases, and their coarsening is wellestablished in other fields such as metallurgy. Phase-field models, especially for alloys, are used in metallurgy to better understand microscopic evolution and changes in the material, as evidenced by recent studies by Geslin et al. (2015), Radhakrishnan et al. (2018), and Ansari et al. (2021). The phase-field method is a numerical simulation technique used to describe phase boundaries and microstructures within materials. Phase-field models are a powerful tool in materials science, physics, medicine, biology, and earth sciences for simulating morphological changes in materials. These simulations use numerical methods to represent the spatial and temporal evolution of individual phases, including their interfacial motions, and can consider different characteristics such as chemical compositions and crystallographic orientations. The microstructure results from the shape and volume of the individual phases, which in turn are formed by diffusive, mechanical, thermal, electrochemical or magnetic driving forces [126]. Phase-field models first gained prominence in the 1970s and 1980s when scientists explored phase transformations, specifically focusing on solidification and microstructure evolution. General background 29 Pioneers like J.W. Cahn and J.E. Hilliard made significant contributions, laying the foundational groundwork for understanding the mathematics of phase separation and transitions. In the 1990s, researchers such as W.J. Rappel and A. Karma expanded these studies to include solidification and dendritic growth. With recent advancements in high-performance computing, phase-field models now enable detailed simulations, leading to technological progress through new mathematical models, algorithms, and a broadening range of applications. A defining characteristic of phase-field models is their ability to depict distinct phases within a system using an order parameter. This parameter generally ranges from 0 to 1, with 0 representing one phase, 1 representing another, and values between 0 and 1 indicating the interface that separates these phases. The interface demonstrates diffusivity, creating a gradual transition across a finite area. This inherent diffusivity enables a realistic representation of phase changes as they occur. Additionally, phase-field models streamline the depiction of complex interfacial motions, such as merging, dissolution, and breakup, by reconstructing the interface through a set of continuous field variables. The basic principle of a phase-field model is depicted in Figure 1.13. In a phase-field model, the driving force is characterised by the free energy functional. Different components of the free energy can be included, each incorporating distinct terms tailored to the specific aspects of the energy balance that need representation in the simulation. In general, the free energy functional can be defined as follows [126]: 𝐹 =   𝑓  + 𝑓  + 𝑓  + 𝑓   𝑑𝑉  , (Eq. 5) where fc is the chemical free energy density, fg is the gradient free energy density, fm is the mechanical free energy density, and fe is the external free energy density. The chemical free energy density, also known as homogeneous or bulk free energy density, promotes phase separation in the absence of an interface, while the gradient energy density penalises the formation of sharp interfaces by imposing penalties. The mechanical free energy density includes contributions from elastic displacements caused by stress and strain. Additionally, all external forces, such as electrostatic or magnetic influences, are incorporated into the external free energy density. During simulations, concerted efforts are made to minimise the system's free energy, which leads to changes in the phases within the system. General background 30 Additionally, phase-field models also have the capacity to describe multiple phases concurrently as multi-phase systems. Figure 1.13: Principle of the order parameter (η) in phase-field modeling. The phases are captured by values between 0 (black) and 1 (white). Modified from [127]. The relevance of phase-field models in the pharmaceutical environment is increasingly evident, with a growing body of literature addressing the microscopic behaviour of substances within this context. Van der Sman (2016) enhanced a two-dimensional phase-field model with theories on thermodynamics and diffusion kinetics to simulate ice crystal growth in a sugar solution during freezing. Similarly, Fan et al. (2018) developed a sophisticated twodimensional phase-field model to describe the phenomenon of freeze-concentration during freezing, based on first-principle assumptions. Li and Fan (2020) modelled macroscopic freezing in a cylindrical vessel, while Li et al. (2022) applied phase-field modeling to explore dendritic morphologies and ice crystal growth inhibition in a sucrose solution. The investigations in this work are based on the model used by Mukherjee et al. (2009), which was originally used by the authors to elucidate the effect of misfit strain and interface curvature on the growth of a single precipitate in a supersaturated matrix. A detailed description of the model can be found in the Results section of this work. Materials and methods 37 analyses were also carried out and are included in Appendix A. The simulation results were visualised using the open-source data analysis software ParaView (Sandia National Laboratories, Kitware Inc, Los Alamos National Laboratory). 2.2.6 Image analysis To determine recrystallisation rates in annealing experiments as well as investigate the freezing behaviour of disaccharide solutions, individual ice crystals on FDM images were marked in the open-source scalable vector graphics editor Inkscape (Inkscape Community). The automatic counting function in Inkscape was then used to determine the number of ice crystals in each image. This process is shown in Figure 2.6. Figure 2.6: Counting of individual ice crystals from an FDM image section of a 10 % (w/w) trehalose solution during annealing at -6 °C for 10 min (left) and 60 min (right). The blue dots were manually placed to mark each individual ice crystal. Results and discussion 38 3 Results and discussion The findings of this thesis are presented in the following sections. As described in the Objectives chapter, the main aspects can be summarised as follows. The first part focuses on the investigation of the microstructure consisting of a particulate (here: ice crystal) and matrix phase during freezing of an aqueous disaccharide solution. In the second part, a phase-field simulation was implemented and calibrated with experimental data to model the coarsening behaviour of the microstructure in the frozen solution during annealing. The final step of this work was to compare and discuss data obtained from the literature with data generated from the phase-field simulations. 3.1 Microstructure formation during freezing In Chapter 1.5, the formation of microstructure during the freezing of an aqueous solution, as often discussed in the literature, was explained. It was also noted that samples develop a complex thermal history during the freezing process, namely, the temperature of the sample does not merely follow the shelf temperature (see Chapter 1.4.1). This raises the question of how much the thermal history, derived from the measured temperature profile during freezing, influences the formation of the microstructure. To assess this, the temperature profiles of samples during freezing in a lyophiliser were analysed, followed by an investigation of the impact of thermal history on the microstructure in a frozen sucrose solution using a freeze-drying microscope. Given that disaccharide solutions ranging from 5 % to 20 % (w/w) are commonly used in the field of lyophilisation [70,83,71], a 10 % (w/w) sucrose solution was employed for this work. 3.1.1 Recalescence and crystallisation in samples during freezing To investigate the thermal history of samples during freezing, vials were filled with a 10 % (w/w) sucrose solution and placed in a lyophiliser. Thermocouples were positioned close to the glass bottom inside the samples to record the temperature. After a short precooling phase to equilibrate the sample temperature to approximately 5 °C, the shelves were cooled to -40 °C at a rate of 1.2°C/min. Images of the filled vials were taken continuously during this process. The results are shown in Figure 3.1. Results and discussion 39 Figure 3.1: Shelf and sample temperature during freezing at 1.2 °C/min. The sample temperature was measured near the bottom of the vial. Pictures were taken of samples at different stages: supercooling (A), nucleation and recalescence (B), crystallisation (C), and complete solidification (D). During the cooling of the shelves, the temperature inside the samples also decreased. Notably, a delay between shelf and sample temperature became apparent due to the imperfect conduction of heat and the placement of the thermocouple with a small distance of approximately 5 mm from the bottom of the vial. The sample temperature decreased to nearly -13 °C without freezing, which means that the solution became supercooled once falling below its equilibrium freezing temperature (A). Upon the formation of stable nuclei in the supercooled solution, the initiation of ice crystal growth lead to a sudden increase in the sample temperature as a result of the heat of crystallisation (B). Once the sample reached its maximum temperature during recalescence (B), the temperature remained elevated for a period of time (B→C). This is due to the heat balance between the latent heat generated by the freezing water in the sample and the heat removed by the cooling shelf of the lyophiliser, Results and discussion 40 as will be explained in the following. The total amount of latent heat generated by water during freezing is approximately 79.7 cal/g [66]. This means that about 79.7 calories of heat are generated to completely freeze one gram of water. In the experiment, however, this heat was not immediately removed from the samples. Instead, the supercooled sample took up only a portion of that heat, which in turn raised its temperature close to its equilibrium freezing temperature (B→C). Depending on the degree of supercooling, approximately 15 cal/g of heat are nearly immediately taken up by the sample as sensible heat [2]. This indicates that only a portion of the water was rapidly frozen once nucleation occurred. Subsequently, the temperature of the sample was continuously reduced via the cooling shelves, which in turn allowed more water to crystallise and generate more latent heat. In order to completely solidify the sample, the totality of latent heat of ice crystallisation (79.7 cal/g) had to be removed via the cooling shelves of the freeze-dryer throughout the entire freezing step once nucleation had occurred (B→C→D). As soon as the water content available for freezing was completely crystallised, and no additional heat could be generated, the sample temperature decreased and converged almost to the shelf temperature (D). In summary, the “freezable” water in the sample did not freeze directly but necessitated the removal of the heat it generated while further crystallising. The time required to achieve this caused a delay in the sample temperature drop relative to the shelf temperature. This delay resulted in a certain time the sample remained above its nucleation temperature (B→C). This entire process was also observed via a camera in the freeze-dryer. As a side note, the placement of the samples for the image recording did not follow the scheme previously shown in Chapter 2.2.1. Instead, the samples were positioned (without any thermal sensors) directly in front of the camera in order to take pictures. The supercooled and transparent sample (A) became slightly opaque once nucleation occurred (B). Over the course of freezing a gradual increase in opaqueness was observed from the bottom to the top of the vial (C), until the sample became completely opaque after complete solidification (D). It can be inferred that the opaqueness of the sample is dependent on the amount of ice in the sample. More specifically, the liquid phase and the crystalline ice phase have different refractive indices, which means that the total surface area of ice crystals, where refraction of light occurs, defines the opaqueness of the sample. During recalescence, the volume fraction of the crystalline phase is lower than in a completely solidified sample due to the increased sample Results and discussion 41 temperature. Therefore, the total surface area of the crystalline phase is also lower than in a completely solidified sample. Furthermore, the gradient of opaqueness during freezing (C) can be explained by the directional removal of heat from the sample (and therefore further crystallisation) at the bottom of the vial by the cooling shelves. This phenomenon is particularly of interest, because it indicates that the temperature conditions in the sample are not uniform but rather dependent on the distance to the cooling shelf. 3.1.2 Characterisation of the thermal profile during freezing The temperature profile shown in Figure 3.1 is typical of lyophilisation samples, and it is often used to describe the phenomenon of supercooling during freezing [2,16,44,86]. Several studies have been carried out to investigate the thermal evolution of samples during lyophilisation, focusing on the influence of the temperature ramp of the shelves and the nucleation temperature of the samples [86,134–136]. These studies typically compare the detected nucleation temperature and freezing rate with primary drying times. The underlying assumption is that the nucleation temperature determines the number of stable nuclei that will grow into ice crystals and later form the pores in the lyophilisate (see Chapter 1.5). This seems reasonable since the amount of water available for freezing is finite, which means that the higher the number of ice crystals, the smaller the average ice crystal size must be. It is also assumed that this pore size will affect the drying time of the sample. Specifically, smaller pore sizes should result in increased resistance to vapour flow through the lyophilisate, thus reducing the rate of sublimation during drying. However, these studies tend to overlook the fact that during lyophilisation the samples are exposed to elevated temperatures due to both recalescence and the subsequent crystallisation phase until the total heat of crystallisation is removed from the sample. For the purposes of this work, the time during which the temperature of the samples is increased as a result of recalescence and crystallisation is referred to as the residence time (tr). The residence time can be defined for different temperature limits, i.e., when does the sample temperature exceed or fall below a certain temperature, as shown in Figure 3.2. In addition, the area under the curve (AUC) can be derived from the temperature profile and the corresponding temperature limit of the residence time. It should be noted that neither the Results and discussion 42 residence time nor the AUC for the respective residence time are physical values intended for any kind of thermodynamic quantification, but serve in this work as a measure to evaluate the effect of freezing on the exposure of the sample to elevated temperatures. Figure 3.2: Temperature profile of a 10 % (w/w) sucrose solution during freezing with nucleation temperature (Tn) and residence times (tr) for temperature limits of -5 °C and -10 °C. The area under the curve (AUC) is also shown for temperature limits of -5 °C (red) and -10 °C (red and blue). To investigate the residence time during freezing, vials were filled with 3 ml of a 10 % (w/w) sucrose solution and placed in the lyophiliser. Two empty vials were left between each sample to reduce heat transfer between samples, as discussed in Chapter 2.2.1. Thermocouples and wireless temperature sensors were placed approximately 1 mm above the bottom of each sample vial to ensure comparability of individual measurements, as a thermal gradient is expected in the samples due to cooling from the bottom of the vials. It should be noted that the use of thermocouples/wireless temperature sensors introduced additional surfaces where nucleation could occur. However, this could not be avoided as the only way to measure temperature in this work was by this invasive method. The samples were pre-cooled to approximately 5 °C and then frozen to -40 °C using two different freezing rates. The first rate was 1.2 °C/min, which is close to the upper limit possible with the freeze-dryer in this work, while the second rate of 0.1 °C/min served as a very slow freezing rate. The nucleation temperature was first determined for each measurement and is shown in Figure 3.3. Results and discussion 43 Figure 3.3: Effect of freezing rate on nucleation temperature (Tn) in a 10 % (w/w) sucrose solution. Standard deviations (error bars) were determined for thermocouples with n = 8 and wireless sensors with n = 16. The nucleation temperature did not show a significant dependence on the freezing rate for both thermocouple and wireless sensor measurements (p > 0.05). This is consistent with the findings of Searles et al. (2001), who reported no effect of freezing rate on nucleation temperature in the range of 0.05 °C/min to 1 °C/min. It can be concluded that the higher freezing rate was not sufficient to lower the temperature fast enough to affect nucleation. Instead, the supercooled samples all became unstable at around -11.9 °C ± 2.1 °C and started to form the first stable nuclei. The residence time was then examined in relation to the freezing ramp and is shown in Figure 3.4. For this purpose, temperature limits between -2 °C and -10 °C were chosen and the corresponding residence times were extracted from the temperature profiles. The lower temperature limit of -10 °C was chosen for the current study because it is close to the average nucleation temperature of the 10 % (w/w) sucrose solution in the freeze-dryer, while the upper limit of -2 °C is only slightly below the equilibrium freezing temperature of the solution. In addition, temperature intervals of 2 °C were chosen to more accurately capture the temperature profile observed during recalescence and crystallisation. It should be noted that this approach is only an attempt to characterise the thermal profile and does not include any thermodynamic considerations regarding the freezing solution. Results and discussion 44 Figure 3.4: Effect of freezing rate on residence times (tr) in a 10 % (w/w) sucrose solution. Standard deviations (error bars) were determined for thermocouples with n = 8 and wireless sensors with n = 16. Significant differences in residence times were observed between freezing ramps of 0.1 °C/min and 1.2 °C/min for each temperature limit (p < 0.05). At the highest freezing rate of 1.2 °C/min, which is close to the maximum capacity of the equipment and therefore represents the upper limit of heat removal possible in a lyophiliser, the residence time was between 6.3 min ± 3.7 min and 9.5 min ± 1.1 min, depending on the temperature limit. The lower freezing rate of 0.1 °C/min resulted in residence times between 13.0 min ± 3.7 min and 23.0 min ± 6.1 min, also dependent on the temperature limit. It can be concluded that the freezing ramp can be used to influence the residence times of the samples during freezing. An alternative method of assessing the effect of freezing on the exposure of the sample to elevated temperatures was by determining the AUC of the recalescence and crystallisation phase using the same temperature limits of -2 °C and -10 °C. This approach considers the nonlinear nature of the temperature profile during freezing, and might therefore capture details better than the residence time. The results of the evaluation are shown in Figure 3.5. Significant differences were observed between the freezing ramps of 0.1 °C/min and 1.2 °C/min (p < 0.05), except for a temperature limit of -2 °C (p > 0.05). The AUC showed similar trends to the evaluation of the residence time. However, the AUC exhibited a linear increase with decreasing temperature limit, whereas the residence time showed a flattened profile. Results and discussion 45 This difference can be attributed to the non-linear temperature profile during recalescence and crystallisation Figure 3.5: Effect of freezing rate on the area under the curve (AUC) of the thermal profile of a 10 % (w/w) sucrose solution during freezing. Standard deviations (error bars) were determined for thermocouples with n = 8 and wireless sensors with n = 16. In the field of lyophilisation, there is a consensus that the nucleation temperature affects the primary drying time [86,137–140]. Specifically, higher nucleation temperatures are associated with shorter primary drying times due to fewer stable nuclei and consequently larger average ice crystal sizes. This observation aligns with classical nucleation theory, which explains the relationship between nucleation temperature and nucleation rate (see Chapter 1.5). However, previous studies have largely overlooked the role of the recalescence and crystallisation phase, focusing instead on attributing ice crystal size prior to drying to nucleation kinetics. The goal of this work is not to dispute classical nucleation theory but to explore whether the temperature increase from the latent heat of crystallisation impacts the microstructure and if its significance has been underestimated. This first raises the question whether nucleation temperature has an impact on the thermal profile of a sample during freezing. To investigate this, the residence times and AUC from the freezing experiments were plotted against the nucleation temperatures for each sample. The results are shown in Figure 3.6. For clarity, only data pertaining to a temperature limit of -6 °C with wireless sensors are displayed. Additional data can be found in Appendix B. Results and discussion 46 Figure 3.6: Comparison of residence time (tr) and area under the curve (AUC) of the recalescence and subsequent crystallisation phase with the nucleation temperature (Tn) during freezing of a 10 % (w/w) sucrose solution measured with wireless temperature sensors. The residence time and the AUC of the recalescence both exhibited a positive correlation with the nucleation temperature. Although the coefficient of determination R² is rather low, it has to be considered that uniform temperature conditions cannot be guaranteed in a freeze-dryer, and that the vials were distributed over nearly the entire shelf (see Chapter 2.2.1). Nevertheless, the quality of the data is comparable to other works in this field [86]. The observed dependence between recalescence and subsequent crystallisation, here quantified by residence time and the AUC, and the stochastic nucleation temperature supports the following hypothesis about the formation of the microstructure in a frozen solution. In the following, a hypothesis will be made that explains the relationship between the thermal profile of a sample during freezing and the nucleation temperature. During the freezing step, samples typically become supercooled as described in Chapter 1.4.1. As soon as nucleation occurs, latent heat of crystallisation is generated during the recalescence and crystallisation phase and must be removed by the cooling shelf. However, the amount of latent heat released immediately after nucleation depends on the degree of supercooling. This phenomenon is illustrated in Figure 3.7. The initial uptake of the latent heat of crystallisation by the sample is only possible until the sample temperature nearly reaches its equilibrium freezing temperature (Tf). Exceeding this temperature would mean that crystallisation could produce enough heat to melt the crystalline phase completely, which would be absurd from a thermodynamic point of view. Results and discussion 53 have a concave bottom, which results in contact with the cooling shelf only at the outer edge of the vial, further reducing the contact area. For the subsequent experiments, 1 µl of a 10 % (w/w) sucrose solution was pipetted onto a glass plate, covered with a second glass plate, and placed into the freeze-drying stage of the microscope (refer to Chapter 2.2.2). The temperature of the cooling block was rapidly reduced to -40 °C at a rate of 50 °C/min to freeze the sample. This steep freezing rate and the low volume-to-contact area ratio were aimed at minimising changes in the microstructure that could be enabled by a rise in the sample temperature due to heat of crystallisation, thereby addressing the issues of temperature gradients and heterogeneity typical of a vial setup. Experiments were performed with and without a spacer between the glass plates to determine the most suitable sample preparation. Once completely frozen, the temperature was raised to -5 °C at a rate of 10 °C/min and maintained to anneal the samples, simulating the temperature conditions during the recalescence and the subsequent crystallisation phase. The results of these experiments are shown in Figure 3.11. Figure 3.11: Freeze-drying microscopy of a 10 % (w/w) sucrose solution with (top) and without (bottom) a spacer after annealing at -5 °C for various annealing times (ta). The spacer increased the sample height to approximately 70 μm. The frozen sucrose solution showed a lack of well-defined structures under the freeze-drying microscope when a spacer was used. Even after a short annealing phase of 5 minutes, during which the particulate ice phase was expected to coarsen due to Ostwald ripening, the Results and discussion 54 identification of individual structures remained speculative. In contrast, omitting the spacer facilitated the formation of a thin layer, with a calculated sample height of approximately 5 µm. This setup allowed for the identification of distinct ice crystals with grain-like shapes. This concept is further illustrated in Figure 3.12. Figure 3.12: Freeze-drying microscope image of a 10 % (w/w) sucrose solution after annealing at -6 °C for 10 min (left) and schematic depiction of the cross section of the sample (right) with bordering ice crystals (A) or ice crystals separated by the matrix phase (B). The clear visibility of individual ice crystals can be attributed to their simultaneous contact with both the lower and upper glass plates. This arrangement causes the phase boundary between the ice crystal and the amorphous matrix to extend across the entire distance from one glass plate to the other. In cases where ice crystals are sufficiently small to allow multiple crystals to fit between the glass plates, phase boundaries would form parallel to the plates, hindering the identification of individual ice crystals. Notably, when the frozen sample was not subjected to an additional annealing step, even without a spacer, the visual characteristics of the solution immediately after freezing did not align directly with the presence of grain-like ice crystals as observed after the annealing process, as shown in more detail in Figure 3.13. Instead, a fine texture was observed, and the determination of its microstructural properties was not possible due to the insufficient magnification of the freeze-drying microscope. The images of non-annealed samples shown in Figure 3.13 prompt an investigation into the arrangement of the crystalline and amorphous components within the frozen sample. Two hypotheses can be considered. The first suggests that, similar to the annealed samples, grainlike ice crystals might be present within a continuous matrix phase, but significantly smaller in Results and discussion 55 size. The overlap of these small ice crystals, in addition to their small size, may hinder their clear observation using a freeze-drying microscope. Alternatively, it is plausible that the ice crystals in non-annealed samples form differently, not appearing as grain-like structures. Alternatively, it is plausible that the ice crystals in non-annealed samples are present in a different form rather than grain-like ice crystals. This phenomenon is intriguing because of the significant size difference between the structures observed in non-annealed samples under the freeze-drying microscope and the pores typically found in lyophilised products. As mentioned in Chapter 1.5, different studies using different methods have reported average pore diameters of approximately 20 µm to 60 µm in non-annealed 5 % (w/w) to 20 % (w/w) sucrose samples [70,83,98,71]. This order of magnitude is comparable to that observed in annealed samples from the freeze-drying microscope (see Figure 3.13). However, this does not hold true for non-annealed samples with their fine structures. Figure 3.13: Freeze-drying microscopy of a 10 % (w/w) sucrose sample without (left) and after 10 minutes of annealing at -6 °C (right). The freezing temperature was approximately -17 °C. The next step involved testing different freezing ramps in the freeze-drying microscope to identify at what temperature nucleation occurred and whether the nucleation temperature could be affected by the freezing ramp in this experimental setup. Nucleation temperatures were determined by capturing images as the sample cooled until crystallisation was visible. This transition was observed in the camera view as texture-like fine structures began to form, as depicted in Figure 3.13. Due to camera limitations at rapid freezing rates, images were captured every 1 °C, resulting in a nucleation temperature data resolution of 1 °C. Figure 3.14 displays the nucleation temperatures for different freezing rates observed with the freezedrying microscope. Results and discussion 56 Figure 3.14: Effect of freezing rate on nucleation temperature (Tn) in a 10 % (w/w) sucrose solution under a freeze-drying microscope. Standard deviations (error bars) were determined from triplicates. The average nucleation temperature in the freeze-drying microscope was determined to be -16.9 °C ± 1.1 °C. No significant differences in nucleation temperatures were observed between freezing ramps from 0.1 °C/min to 20 °C/min (p > 0.05). In contrast to the freezing ramps used in the lyophiliser, steeper ramps were tested in the freeze-drying microscope to determine whether the highest possible freezing ramp in a lyophiliser of 1.2 °C/min was simply insufficient to affect the nucleation temperature. Consequently, the freezing ramp was found to be ineffective in influencing the nucleation temperature in both lyophiliser and freezedrying microscopy experiments. However, the nucleation temperature in the lyophiliser was found to be -11.9 °C ± 2.1 °C, which significantly differed from the values in the freeze-drying microscope (p < 0.05). To ensure comparability between freeze-drying microscopy and lyophilisation experiments, the nucleation temperature in the freeze-drying microscope had to be adjusted accordingly. Efforts were made to induce nucleation at approximately -12 °C, including attempts to hold the temperature for a prolonged time, abruptly reducing and increasing the pressure in the freeze-drying stage, and oscillating temperature within ± 1 °C of the desired nucleation temperature. The most reproducible results were obtained when nucleation was induced by repeated gentle taps on the sample holder while maintaining the sample temperature at -12 °C. Further examination of the samples was carried out in the freeze-drying microscope using polarised light, as shown in Figure 3.15 for frozen and non-annealed samples. The Results and discussion 57 principles of polarisation and birefringence have already been described in Chapter 2.2.3. The samples were cooled to -12°C and nucleation was induced by tapping the sample holder. The samples were then cooled to -40°C at 10°C/min and images were taken. Figure 3.15: Freeze-drying microscopy of a 10 % (w/w) sucrose sample with (right) and without (left) polarisation at -40 °C with preceding nucleation at around -12 °C. Results and discussion 58 The freeze-drying microscope images show identical sections with and without polarised light. In addition, the outer edge of the frozen sample was deliberately included in the 100x and 200x magnifications (bottom right-hand corner of each image). Due to the high freezing rate of 10 °C/min, the samples were cooled to -40°C within 3 minutes of nucleation. No separate annealing step was performed for these samples. The polarised light microscopy images also indicate the presence of feather-like structures without a clear pattern similar to an ice flower on a window. Additionally, the structures seem to have originated and spread out from the phase boundary of the solution. To investigate the transition from ice crystals in non-annealed samples to grain-like ice crystals, the sample was heated from -40 °C to -6 °C at 10 °C/min and annealed for 120 minutes. The corresponding images are shown in Figure 3.16. Irregular structures were observed in the non-annealed sample at the top of Figure 3.16, while single grain-like ice crystals became visible after short annealing times (ta = 10 min). Each individual ice crystal appeared to have a single uniform colour. It is also noticeable that neighbouring ice crystals often exhibited the same colouration. These clusters of ice crystals of the same colour were often located where irregular structures of a similar colour had previously been present. In order to further investigate the effect of the recalescence and crystallisation phase on the microstructure in a frozen solution, the annealing step in the freeze-drying microscope was designed to mimic the conditions experienced by the samples in the lyophilisation experiments from Chapter 3.1.2. However, the temperature profile of the sample in the lyophiliser during freezing was not linear or stepwise, so it would not have been feasible to implement such a profile in the freezedrying microscope. In addition, care must be taken to ensure that the sample temperature does not exceed the equilibrium freezing temperature of the solution, as this would result in complete thawing of the sample. Some temperature fluctuations are to be expected in the freeze-drying microscope and therefore it would have been too risky to use a set temperature slightly below the equilibrium freezing temperature for a 10 % (w/w) sucrose solution of about -0.6 °C [66]. Therefore, in the following experiments, an annealing temperature of -6 °C was chosen, which, based on previous experience, clearly demonstrated the effect of annealing on the microstructure, but was sufficiently far from the equilibrium freezing temperature to prevent the complete thawing of the sample. The residence times in the lyophilisation experiments from Chapter 3.1.2 were taken into account for the annealing time for the freeze- Results and discussion 59 drying microscopy experiments. For a temperature limit of -6 °C, the residence times for freezing ramps of 1.2 °C/min and 0.1 °C/min were rounded off to 10 minutes and 20 minutes, respectively. Figure 3.16: Freeze-drying microscopy with polarised light of a 10 % (w/w) sucrose sample during annealing at -6 °C. Results and discussion 60 The microstructure of a frozen solution was assessed by counting individual ice crystals after mimicking the temperature conditions of the recalescence and crystallisation phases experienced during freezing in a lyophiliser. As previously discussed in Chapter 1.5, the number of ice crystals serves as an indicator of the eventual pore size after drying. This is because the freezable water in the solution is distributed among all the ice crystals formed during freeze-concentration, i.e., a larger number of ice crystals corresponds to smaller average ice crystal sizes. For this evaluation, samples were either cooled using different freezing ramps until nucleation occurred or nucleation was induced at -12 °C. Following nucleation, the samples were rapidly cooled to -40 °C at a rate of 50 °C/min to minimise structural changes that might occur above the glass transition temperature. The samples were then heated to an annealing temperature of -6 °C at 10 °C/min and held for 20 minutes. Images were taken after 10 minutes and 20 minutes, corresponding to residence times of freezing ramps of 1.2 °C/min and 0.1 °C/min in a freeze-dryer. The number of ice crystals was quantified through image analysis, as detailed in Chapter 2.2.6, across a viewing area of approximately 1.12 × 10³ μm². The results are shown in Figure 3.17. Figure 3.17: Effect of freezing rate and annealing time (ta) on the number of ice crystals in a 10 % (w/w) sucrose solution. Standard deviations (error bars) were determined from triplicates. For all freezing ramps, as well as for enforced nucleation at -12 °C, no significant difference was observed in the number of individual ice crystals (p > 0.05), indicating that under uniform Results and discussion 61 annealing conditions, the coarsening of the crystalline phase occurred independently of the preceding freezing conditions. However, a significantly lower number of ice crystals was counted after 20 minutes of annealing compared to 10 minutes, across all freezing ramps and enforced nucleation conditions (p < 0.05). The relationship between the number of ice crystals and the nucleation temperature is illustrated in Figure 3.18. It is important to note that the nucleation temperature for each measurement was determined with an accuracy of 1 °C, due to limitations in capturing more images with the camera during the steep freezing ramps. Figure 3.18: Effect of nucleation temperature (Tn) and annealing time (ta) on the number of ice crystals in a 10 % (w/w) sucrose solution. These results also show that over a range of nucleation temperatures from -12 °C to -18 °C, the resulting number of individual ice crystals did not change significantly after the samples were exposed to either 10 minutes or 20 minutes of annealing (p > 0.05). However, there was a significant difference between the datasets corresponding to the two different annealing times at any given nucleation temperature (p < 0.05), similar to the results from Figure 3.17. The results of the annealing experiments in the freeze-drying microscope can be summarised as follows:  Freezing a small volume of sample while preventing the temperature from rising due to recalescence and crystallisation leads to the formation of fine feather-like structures. Results and discussion 62  If these structures are composed of small individual ice crystals, then these ice crystals are one or more orders of magnitude smaller than any pore found in a lyophilisate.  When these frozen samples are annealed in a freeze-drying microscope under conditions that mimic typical residence times in a lyophiliser, individual grain-like ice crystals with sizes comparable to pores in a dry lyophilisate are formed.  The quantity (and therefore size) of these ice crystals depends mainly on the conditions during the recalescence and crystallisation phases, and not on the nucleation temperature. It is important to note that the reduction in sample thickness can impact the rate of Ostwald ripening of ice crystals during annealing. The spatial confinement between glass plates alters the shape of ice crystals by flattening them, as illustrated previously in Figure 3.12. Ostwald ripening, described in Chapter 1.4.2, is a curvature-dependent process and may decelerate when fewer curved surfaces are present. Therefore, it is crucial to clarify that the objective of these experiments is not to determine an exact pore size based on freezing and annealing conditions. Rather, the goal is to assess the qualitative effect of recalescence and the subsequent crystallisation phase on the morphological changes within the microstructure of the frozen solution. 3.1.5 Impact of thermal history during freezing on the microstructure The following hypothesis is proposed to interpret the experimental data from previous chapters, with additional context on the freezing of a glass-forming disaccharide solution. The growth of ice crystals from a nucleus in a solution can be classified into three main types: dendrites, irregular dendrites, and spherulites, influenced by factors such as freezing rate, solute composition, and concentration [48]. Regular dendrites form when there is sufficient time for ice crystal growth, while faster crystallisation leads to the development of irregular dendrites and spherulites [2]. Given that the samples experienced supercooling both in the freeze-drying microscope and in the lyophiliser, it is likely that the crystallisation process occurred relatively quickly after nucleation. This rapid formation is evident in Figure 3.15, where the structures appear less organized than, for instance, a snowflake, which exemplifies dendritic ice crystal growth. Structures like dendrites, irregular dendrites, and spherulites are Results and discussion 69 setup of a freeze-drying microscope. The concept of planes of curvature is illustrated in Figure 3.22. Figure 3.22: Schematic representation of the morphology of ice crystals between two spatial boundaries. If the distance between the boundaries is sufficiently large, spherical ice crystals with two planes of curvature are formed (A), whereas if the distance is small, cylindrical ice crystals with only one plane of curvature are present (B). Ice crystals in three-dimensional space without any boundaries exhibit two curvature planes, namely a horizontal and a vertical plane. However, when the distance between the boundary surfaces that limit the ice crystals becomes smaller than the diameter of the ice crystals, then the ice crystals start to deform until they assume cylindrical shapes that only have a single curvature plane. The driving force behind Ostwald ripening, the phenomenon that is of interest for this work, is based on the Gibbs-Thomson effect, which refers to the dependence of vapour pressure or chemical potential on the interfacial energy/curvature of a surface [146]. As explained in Chapter 1.4.2, the saturation vapour pressure deviates from a purely temperature-dependent material property when the material is present in sufficiently small particles such as microscopic ice crystals. A mathematical relationship can be established between curvature and the deviation in saturation vapour pressure, which is described in Appendix A for phasefield models. If a circle (two-dimensional) and a sphere (three-dimensional) with radius r are considered, the mean curvature is 1/r and 2/r, respectively. This means that the deviation in saturation vapour pressure doubles when a circle is expanded into a sphere, i.e., another plane of curvature is added (see Equation App. 1). Accordingly, phenomena such as Ostwald ripening would also proceed at a different rate for a circle and a sphere of the same radius. Results and discussion 70 In freeze-drying microscopy experiments, the sample volume is confined to a thin film between two glass plates. The distance between the two glass plates can be controlled using a spacer. If no spacer is used, the distance between the two glass plates will adjust according to the sample volume (see Chapter 3.1.4). When the distance between the glass plates is sufficiently small compared to the size of the ice crystals, the ice crystals are forced into a cylindrical shape (see Figure 3.22). Subsequently, the plane of curvature in the vertical direction becomes negligible, reducing the system to a single plane of curvature. In this state, annealing still causes Ostwald ripening to occur, but the curvature-dependent coarsening of the ice crystals is expected to be altered due to the reduction of the system to a single plane of curvature. Precisely this property was utilised to carry out a simulation-based investigation of the coarsening, because two-dimensional simulations are representative of the freeze-drying experiments due to the single curvature plane, while three-dimensional simulations with their two curvature planes reflect the ice crystals in a lyo vial. Phase-field modelling was chosen as the simulation type for the following reasons. Phase-field models can be easily extended from two-dimensional to three-dimensional, while the overall mathematical basis and simulation parameters remain identical. In addition, phase-field models inherently incorporate the concepts of curvature (interfacial energy) and saturation vapour pressure (composition/chemical potential). They provide a straightforward way of describing the temporal evolution of phases in such a system and are often used to predict coarsening phenomena. The simulation strategy is illustrated in Figure 3.23 and explained below. First, samples for freeze-drying microscopy experiments were prepared with small volumes in order to reduce the sample height as much as possible. The aim was to reduce the vertical plane of curvature to such an extent that its contribution to the coarsening of the ice crystals became negligible. In Chapter 3.2.2.2.3, annealing experiments were carried out with these samples in the freeze-drying microscope to determine the recrystallisation rates for a single plane of curvature at different annealing temperatures. This experimental data was used to calibrate a phenomenological two-dimensional phase-field model in order to determine simulation parameters that would capture the coarsening behaviour as accurately as possible (see Chapter 3.2.4). Results and discussion 71 Figure 3.23: Approach to the prediction of morphological properties of the microstructure in a lyophilisate in this work: sample preparation for freeze-drying microscopy (1), calibration of a two-dimensional phase-field model (2), determination of simulation parameters (3), prediction of microscopic morphological parameters (4). A three-dimensional phase-field simulation was then implemented using the simulation parameters determined from the two-dimensional simulation in Chapter 3.2.5. The aim of this step was to re-introduce the horizontal plane of curvature, thus reversing the reduction in the number of planes of curvature from the experimental setup. This would then allow Results and discussion 72 conclusions to be drawn for the impact of the annealing process on a frozen disaccharide solution in a lyo vial. Finally, with the three-dimensional extension, morphological parameters were predicted in their temporal evolution, such as the ice crystal size distribution, the surface area of the matrix phase and the conjunction area between ice crystals. The details of each morphological parameter are given in the relevant chapters. 3.2.2 Implementation and calibration of the phase-field model In the following sections, a phase-field model was implemented as a numerical method to investigate the microstructural changes in a frozen solution during annealing, and experimental data were collected to calibrate the relevant model parameters. 3.2.2.1 Formulism of the multiphase-field model The considered phase-field model is based on a particulate phase p (here: ice crystals) and a continuous matrix phase m. These phases are described in space r and time t by the order parameter η(r,t) (crystalline ice phase = 1, amorphous matrix phase = 0) and their composition c(r,t) (pure water = 1, pure solute = 0). The transition between the phases is represented in the phase-field model by a smooth interface with values between the corresponding minima and maxima. The evolution of the microstructure in this phase-field model is governed by the Cahn-Hilliard diffusion equation for conserved field variables: 𝜕𝑐 𝜕𝑡 = ∇ 𝑀 ∇ 𝛿𝐹 𝛿𝑐 , (Eq. 6) and the Allen-Cahn relaxation equation for non-conserved field variables: 𝜕𝜂 𝜕𝑡 = − 𝐿 ∇ 𝛿𝐹 𝛿𝜂 , (Eq. 7) where M is the mobility coefficient, L is the relaxation coefficient and F is the free energy of the system [126]. Since Ostwald ripening is a diffusive process, particle-particle interactions such as coalescence must be inhibited. In nature, fusion of ice crystals is prevented through a mismatch of the crystallographic orientations of neighbouring crystals [103]. Results and discussion 73 This concept can also be adopted in a phase-field model, where η(r,t) is extended to the multiphase-field ηi(r,t) and is given by: 𝜕 𝜂  𝜕𝑡 = − 𝐿 ∇ 𝛿𝐹 𝛿 𝜂  ( 𝑟 , 𝑡 ) , (Eq. 8) with i = 1, 2, …, n and the particulate phase is divided between the n elements of the phase fields. The free energy of the system F is composed of the bulk free energy and the contributions from the gradient free energies: 𝐹 ( 𝑓 , 𝑐 , 𝜂  ) =   𝑓 ( 𝑐 , 𝜂  ) + 𝜅  ( ∇ 𝑐 )  + 𝜅  ( ∇ 𝜂  )  𝑑 Ω   , (Eq. 9) where f(c,ηi) is the bulk free energy density, and κc and κη are the gradient energy coefficients for composition and order parameter, respectively. The bulk free energy density f(c,ηi) promotes phase separation whenever gradient and interfacial energies are not present and is given by: 𝑓 ( 𝑐 , 𝜂  ) = 𝑓  ( 𝑐 )  1 −  𝑊 ( 𝜂  )    + 𝑓  ( 𝑐 )  𝑊 ( 𝜂  )   +𝑃𝜂(1−𝜂)     +𝑄𝜂𝜂         , (Eq. 10) where P is a coefficient for the magnitude of the energy barrier between the mand the pphase, Q is a coefficient for the magnitude of the energy barrier between the n particulate phases, and fm and fp are the free energies of the mand the p-phase, respectively. The interpolation function W(ηi) is used in Eq. (7) and is defined as [147]: 𝑊 ( 𝜂  ) =  0 ; 𝜂   ( 10 − 15 𝜂  + 6 𝜂   ) ; 1 ; for 𝜂  < 0 , for 0 ≤ 𝜂  ≤ 1 , for 𝜂  > 1 . (Eq. 11) The free energies fm and fp are given by: 𝑓  ( 𝑐 ) = 𝐴 ( 𝑐 − 𝑐   )  , (Eq. 12) 𝑓  ( 𝑐 ) = 𝐵  𝑐   − 𝑐   , (Eq. 13) where A and B are positive constants, and 𝑐  and 𝑐  are equilibrium compositions for the mand the p-phase, respectively. Results and discussion 74 In this work, the recrystallisation rate in the model was calibrated using the simulation parameters M and L. Both values remained constant, i.e., when M was changed, L was adjusted to the same value, effectively only impacting the recrystallisation rate and not the morphological properties of the particles. If dissimilar values of M and L had been used in a simulation, the ratio of diffusion to relaxation of the particles would have changed, potentially affecting the morphology adopted by the particles during coarsening. The generation of the two-dimensional and three-dimensional initial conditions is discussed later in Chapter 3.2.3. The dimensionless simulation parameters for the phase-field model are listed in Table 3.1. Table 3.1: Dimensionless parameters used in this work. Parameter Value κc 1 κη 1 A 1 B 1 P 1 Q 2 A semi-implicit Fourier spectral scheme was implemented according to Chen and Shen (1998) by treating linear and second-order operators implicitly and residual terms explicitly. Some stability tests and parameter determination tests are provided in Appendix A. 3.2.2.2 Experimental determination of relevant calibration parameters The phenomenological phase-field model employed to predict the effect of annealing on the microstructure of disaccharide solutions, both during annealing and in their maximally freezeconcentrated state, required calibration of relevant simulation parameters. Initially, experimental values were obtained for the volume fractions of ice crystals and the amorphous matrix phase, which vary with sample temperature and thus differ based on whether annealing is performed or if the solution reaches its maximally freeze-concentrated state. Recrystallisation rates during annealing, which are also temperature-dependent, were then determined experimentally. These experiments included both sucrose and trehalose solutions to also consider a second commonly used excipient in freeze-drying. Results and discussion 75 3.2.2.2.1 Mass fractions during annealing As described in Chapter 1.4.1.2, the mass fractions of solute and solvent in the freezeconcentrated matrix phase of a frozen solution depend solely on the sample temperature. This is important because the temperature is increased during annealing above the onset temperature of melting Tm’ at which water molecules dissolve from the crystalline phase into the matrix phase. Consequently, the ratio of solute and solvent in the amorphous phase changes. These mass fractions in the matrix phase follow the liquidus line of the solute-solvent system. For example, if a 30 % (w/w) sucrose solution becomes completely liquid at -1.6 °C, a freeze-concentrated matrix at this exact temperature (regardless of the initial concentration in the solution) will also contain a 30 % (w/w) ratio of sucrose to water. Therefore, the freezing point depression for different solute concentrations can be used to determine the temperature-dependent mass fractions in the freeze-concentrated matrix phase [60]. In the following, the melting temperature Tm of sucrose and trehalose solutions with concentrations between 30 % (w/w) and 55 % (w/w) was determined using differential scanning calorimetry (DSC) with the protocol from Chapter 2.2.4. An example heat flow is shown in Figure 3.24. Figure 3.24: Exemplary thermograms of sucrose solutions with different solute concentrations during melting. The melting peak was used to determine the melting temperature (Tm) of the respective solution. Results and discussion 76 The glass transition during heating can be observed at approximately -33 °C by the shift in the heat flow from one plateau to another. According to Chapter 1.4.1.2, the glass transition temperature of a maximally freeze-concentrated solution (Tg’) should theoretically be identical for all samples, as it is presumed independent of solute concentration. However, observed differences between the solutions can be attributed to the dissimilar thermal histories arising from varying solute concentrations. Given that nucleation temperature is influenced by solute concentration [46], and considering the demonstrated correlation between nucleation temperature and residence time due to recalescence and crystallisation (as shown in Chapter 3.1.2), it can be inferred that the samples underwent different conditions during freezing. After the glass transition, all curves depicted in Figure 3.24 exhibited similar behaviour during heating, although with variations in peak positions and peak areas. These differences in peak area can be attributed to the varying amounts of water present in each solution. A sample with a lower solute concentration contains more water, thereby providing more ice for melting. Since melting is an endothermic phase transition, it appears as a positive heat flow in the thermogram. As mentioned in Chapter 1.4.2, care must be taken when using terms such as melting in the context of solutions. Melting is used here to describe the dissolution of precipitated ice into the matrix phase, which contains sucrose and water, and not the distinct phase change that occurs at a particular temperature when an entire substance changes phase. Specifically, the melting temperature Tm of the solution has been defined in this work as the temperature at which the heat flux peaks. The temperature of the melting peak, indicating the temperature at which the largest fraction of ice in the sample converts to water, is observed to decrease with increasing solute concentration. This relationship is critical for the calibration of the phase-field model and is detailed in Figure 3.25 for both solutes. As the solute concentration increased, the melting temperature for the corresponding sucrose and trehalose solutions decreased. This is consistent with the phenomenon of freezing point depression, where the freezing temperature of the solution is negatively correlated with the amount of solute added to a solvent [149]. Furthermore, the recorded data are in agreement with values reported in the literature [72,150], with no significant difference between the liquidus lines for sucrose and trehalose solutions (p > 0.05). Results and discussion 77 Figure 3.25: Liquidus line of a sucrose-water (left) and trehalose-water (right) system. The melting temperature of 30 % (w/w) to 55 % (w/w) solutions were measured via DSC. Standard deviations (error bars) were calculated from triplicates. A second-degree polynomial trend line was fitted to the data and can be represented as: 𝑦 = − 0 . 0062 𝑥  + 0 . 2188 𝑥 − 2 . 4362 , (Eq. 14) where y is the melting temperature, and x is the solute concentration, with R² = 0.9985. Equation 14 relates the solute mass fraction in the matrix phase to the ambient temperature for sucrose and trehalose solutions with concentrations ranging from 30 % (w/w) to 55 % (w/w). This experimental range corresponds to melting temperatures of -1.6 °C ± 0.0 °C and -9.3 °C ± 0.1 °C, respectively. Thus, Equation 14 can be used to determine the composition of the matrix phase at annealing temperatures within this range. Next, appropriate annealing temperatures were selected to perform coarsening experiments for the calibration of the phase-field model. Annealing should always be performed above the glass transition temperature Tg’ of the solution, otherwise changes in the microstructure will be in the order of mm/year due to the high viscosity of the material [60]. The rate of coarsening, or recrystallisation rate, of the ice phase is dependent on the annealing temperature, with higher temperatures resulting in faster coarsening [72]. However, while higher annealing temperatures shorten the process time, it is important to avoid setting the temperature too close to the melting temperature of the solution as even small temperature fluctuations in the equipment can lead to the complete thawing of the sample. As discussed Results and discussion 78 in Chapter 3.1.4, the melting temperature of a 10 % (w/w) sucrose solution is at approximately -0.6 °C, which, based on experience, allows annealing temperatures of around -2 °C in the freeze-drying microscope without risking the thawing of the sample. Commonly, -5 °C is used in literature to study the effects of annealing on lyophilisate microstructure [70,71], providing a substantial buffer from the melting point while allowing short annealing durations to noticeably impact the microstructure. Therefore, annealing temperatures of -4 °C and -6 °C were chosen for further experiments and simulations, considering the process duration and at the same time maintaining realistic conditions comparable to published studies. In addition, two different annealing temperatures were chosen to demonstrate the temperature dependency of the process. Since -4 °C and -6 °C were within the range of the experimentally determined polynomial from Equation 14, and sucrose and trehalose solutions did not differ significantly in their freezing point depression, the solute to water mass fractions for both disaccharides at -4 °C and -6 °C were calculated to be 0.398 g/g and 0.461 g/g, respectively. 3.2.2.2.2 Volume fractions during annealing/maximum freeze-concentration The next step involved converting the mass fractions into volume fractions for both the ice and matrix phases under annealing conditions and at maximum freeze-concentration. Due to the unavailability of direct measurements of the matrix phase density at -4 °C, -6 °C, and maximum freeze-concentration, literature data for the density of concentrated sucrose and trehalose solutions at 20 °C were utilised. Potential expansion or contraction of the matrix phase due to low temperatures was not considered. In the first-principles simulation work by Fan et al. (2018), the impact of temperature on the density change of the amorphous phase was limited to the contraction of the residual water content, which is particularly small compared to the solute mass fraction at maximum freeze-concentration. Therefore, any deviation in the volume fraction was deemed negligible in this context. The ice phase did not present this issue, as density data for ice under the specified freezing conditions in a lyophiliser were available (see Table 3.3). First, the water mass fractions in the maximum freezeconcentrated matrix phase for both sucrose and trehalose solutions were gathered from various literature sources and are summarised in Table 3.2. Results and discussion 85 surprising, as a difference in recrystallisation rate was previously found in a separate publication [68]. However, it should be noted that the aim of this preceding publication was not to determine an exact recrystallisation rate, but rather to assess Ostwald ripening and glassy state relaxation qualitatively. Accordingly, the annealing experiment was only performed and evaluated once for both solutes, and no standard deviation could be taken. Klinmalai et al. (2017) also found a slight difference between the two solutes during coarsening. However, the quality of the data is questionable as the standard deviation varies greatly depending on the series of measurements. In addition, it can also be observed from the microscopy images that non-cylindrical 3-dimensional crystals form in some experiments. This is particularly problematic as the radius of an ice crystal can no longer be reliably determined and the measurement now strongly depends on the layer thickness of the sample. An increase in the annealing temperature resulted in a significant increase in the recrystallisation rate from an average for sucrose and trehalose of 26.13 µm³/min at -6 °C to 56.80 µm³/min at -4 °C (p < 0.05). The linear relationship observed between the cube of the circle equivalent radius and annealing time suggests that the coarsening behaviour of sucrose and trehalose solutions in the freeze-drying microscope follows the growth law of the LSW theory. However, this observation only validates Equation 3, as not all the conditions of the LSW theory are met within the experimental setup of the freeze-drying microscope. Therefore, the slopes derived from the linear fits in Figure 3.28 were used to directly determine the recrystallisation rate using Equation 3 instead of the first-principles approach from Equation 4. However, this is perfectly adequate as the phase-field simulation in this work is intended to be a phenomenological model anyway, i.e., the simulation should only quantitatively describe the respective phenomenon from the experiment and does not have to be derived from physical laws. All fitting parameters are summarised in Table 3.4. The linear fit of the cube of the circle equivalent radius of ice crystals during the annealing experiment serves as the basis for calibrating the two-dimensional phase-field simulation. Although other fits were also determined, the use of the cube of the circle equivalent radius is a convenient choice as it is already established in the literature and allows for easy comparison between simulation and experiment due to its linear relationship. Results and discussion 86 Table 3.4: Fitting parameter of number of particles (NoP, power law fit), circle equivalent radius (CER, logarithmic fit) and the cube of the circle equivalent radius (CER³, linear fit) at different annealing temperatures (Ta). Solution Ta [°C] NoP CER CER³ a [-] b [1/min] R² a [µm] b [µm] R² k [µm³/min] R² 10 % (w/w) sucrose -4 3895 -0.592 0.9958 4.943 -2.38 0.9919 54.74 0.9849 -6 8815 -0.628 0.9997 4.024 -3.27 0.9695 27.74 0.9973 10 % (w/w) trehalose -4 8482 -0.689 0.9940 5.433 -4.90 0.9930 58.86 0.9914 -6 5970 -0.563 0.9994 3.685 -1.37 0.9862 24.52 0.9891 3.2.3 Initial conditions for two-dimensional and three-dimensional simulations The specification of initial conditions is a crucial step for in silico methods, where relevant parameters must be defined as the starting point of the simulations. For the phase-field model used in this work, the initial placement and size distribution of the particles (=ice crystals) are essential to simulate the subsequent coarsening process. This section describes the method used to generate the distribution and positioning framework of the particulate phase for all two-dimensional and three-dimensional phase-field simulations. Initially, a basic structure for the positioning of individual particles was established. For twodimensional simulations, two simple geometric arrangements were considered: a square and a triangular lattice, as depicted in Figure 3.29. These configurations represent examples of close-packing that maximize the number of equal-sized particles in a given distribution. It is important to note that the choice of the basic geometric structure is not influenced by physical considerations but is instead chosen solely to establish a starting point for creating the initial conditions. The choice of close-packing structure affects the dimensions and particle area fractions. The triangular lattice results in a y-coordinate length of 𝑁𝑦= √3/2∙𝑁𝑥, due to the positioning of the particle centres. This structure also achieves a higher area fraction of the particulate phase, as it utilises the spaces between particles more efficiently. The area fraction for a square close-packed system is approximately 0.79, while for a triangular closepacked system it is about 0.91 [155]. Results and discussion 87 Figure 3.29: Illustration of the close-packing of circles in square (A) and triangular (B) lattices. Both arrangements exhibit an equidistant distribution regarding the centre of the particles. To reiterate, volume and area fractions are used interchangeably for two-dimensional simulations, as they are assumed to be identical in the experimental setup due to the low layer thickness and cylindrical particle shape. This assumption would be invalid if the particles had too much curvature in the vertical height axis. The volume fractions of the particulate phase during annealing were calculated in Chapter 3.2.3.2 and amounted to an average of 0.793 and 0.828 for annealing temperatures of -4 °C and -6 °C, respectively. Therefore, the triangular lattice was chosen for the two-dimensional simulations in this work, as it is easier to reduce the area fraction of a close packing by simply reducing the particle size than by increasing the area fraction of the square lattice. The exact setting of the area fraction in two-dimensional simulations and the volume fraction in three-dimensional simulations will be discussed in the next chapter. For the initial conditions of the three-dimensional simulations, the basic structure for particle distribution needed to be established as well. When extended into three dimensions, the twodimensional triangular structure transitions into a tetrahedral distribution of particles. Two different close-packings based on a tetrahedral building block are possible, as depicted in Figure 3.30. Both close-packings in three-dimensional have similar packing fractions of approximately 0.74 [156]. Consequently, the choice is therefore not based on physical considerations or volume fractions, but rather on ease of implementation. The hexagonal close-packing has a more repetitive positioning and is therefore slightly easier to implement. Analogous to the reduction of the y-coordinate in the triangular two-dimensional case, the z- Results and discussion 88 coordinate also changes in the hexagonal three-dimensional case to 𝑁𝑧=√6/3∙𝑁𝑥. In summary, the particle positioning was chosen to be triangular for two-dimensional simulations and hexagonal for three-dimensional simulations to simplify the setup and implementation process. Figure 3.30: Arrangement of hexagonal and cubic close-packings in different layers. In the hexagonal system, the particles in the third layer are positioned exactly the same as in the first (blue), while in the cubic system the position is slightly shifted (green). Next, the grid size of the simulation domain and the initial number of particles were established based on computational feasibility. The scaling of the total number of grid points in such a system is depicted in Figure 3.31. Figure 3.31: Scaling of the simulation grid size (Nx · Ny · Nz) with the hexagonal positioning of particles. The tetrahedral basic structure results in 𝑁𝑦= √3/2∙𝑁𝑥 and 𝑁𝑧= √6/3∙𝑁𝑥. Results and discussion 89 When determining the domain size, it was crucial to manage the number of grid points to prevent unreasonably long computation times. This consideration is especially important for three-dimensional simulations, where the number of grid points substantially increases due to the addition of the z-coordinate. Consequently, the total size of the three-dimensional domain was used to define the simulation domain. The base grid size was chosen based on experience with an upper limit of approximately 10 million grid points (gp). According to Figure 3.31, a base grid size of Nx = 240 gp was chosen, resulting in a total of just under 9.8 million grid points, with Ny = 208 gp and Nz = 196 gp. The next stage of the process involved the selection of the number of particles to be used for the starting point of the simulation. In the case of two-dimensional simulations, a total of n particles were placed in each coordinate, resulting in a total of n² particles. In the threedimensional setup with hexagonal positioning, a total of 4n³ particles were used. As with the grid size, an issue may arise when the number of particles becomes too large due to the need to increase the number of order parameters in the multiphase-field approach. To reiterate, the order parameter in the multiphase-field model is the equivalent to the crystallographic orientation of ice crystals. This approach allows a defined number of particles to be computed in a single iteration, provided that the particles in the same order parameter are not in close proximity to each other. If particles of the same order parameter come into contact with each other, unintentional interactions such as fusion of neighbouring particles will occur. Therefore, it is important that the initial particles are distributed over a sufficient number of order parameters in order to prevent contact-to-contact interactions between particles within the same order parameter. The advantage of placing multiple particles into the same order parameter is that this reduces the computational time by a factor equal to the average number of particles per order parameter. For example, placing 10 particles per order parameter reduces the computational time by approximately 10-fold. Therefore, the choice of how many particles are placed in each order parameter is a balance between preventing unwanted interactions while decreasing computational effort. Based on experience, n was chosen to be 6 for three-dimensional simulations, resulting in 864 initial particles. All particles were distributed over 86 order parameters, with an average of 10 particles per order parameter. This means that a maximum of 86 iterations had to be calculated for a single time step instead of 864. The particles were generated with a large Results and discussion 90 spacing within the same order parameter to avoid unwanted interactions after slight coarsening. In the case of two-dimensional simulations, n was increased to 10 to have 100 initial particles in order to observe coarsening over a larger number of particles. Due to the short computational times of two-dimensional simulations, 100 order parameters were used, i.e., each particle was placed in its own order parameter. To initiate coarsening in the phase-field simulation, irregularities in the initial particle size, shape, or position had to be introduced, as will be explained below. Figure 3.32 illustrates the process of constructing the initial condition for the twoand three-dimensional simulations. Figure 3.32: Generation of the initial conditions with a triangular base structure for twodimensional simulations (upper row) and a hexagonal base structure for three-dimensional simulations (lower row). A system with equidistant particles of similar size would not coarsen because the interfacial energy or curvature of all particle surfaces would be identical, thus lacking a driving force for change. Introducing irregularities leads to deviations from ideal packing, affecting the expected area and volume fractions. These deviations were, however, are of no consequence Results and discussion 91 as the natural packing fractions of triangular and hexagonal close-packings differed from the calculated fractions. This means that a separate adjustment of the area and volume fractions regardless of initial irregularities needed to be performed anyway. To introduce irregularity, the particle size was first changed from the densest packing, where adjacent particles are in contact, to a size distribution. A slight randomisation of particle positioning was then introduced, moderately varying the distances between particles while maintaining the basic triangular or hexagonal positioning. The choice of Nx, Ny and Nz is based on the computational effort, as the grid size increases exponentially as shown in Figure 3.31. On the other hand, the distances between two grid points dx, dy and dz in the phase field simulations can be chosen more flexibly as they do not affect the computational effort. For the sake of simplicity, only dx will be referred to as the distance between two grid points, as the same value is usually used for dx, dy and dz in a simulation. Increasing dx will result in a larger physical domain of the simulation as the actual edge length of the simulation is calculated using Nx · dx. However, setting dx to a high value may result in interfaces being covered by an insufficient number of grid points. This might lead to grid effects, where the angle of an interface affects its interfacial energy, causing numerical problems during the simulation. An example of this is shown in Figure 3.33. An increase in the distance between grid points dx was accompanied by deformation of the particles when they were tilted (dx = 1.2) or even loss of curvature when dx became too large (dx = 1.5). As dx increased, so did the physical simulation domain, but the interfaces also became narrower and at a critical point it was no longer possible to accurately represent inclinations at the interfaces. The choice of dx was therefore a trade-off, as lower dx required higher Nx to capture the same physical space, while higher dx reduced accuracy and even led to errors. After careful consideration, a value of dx = 0.9 was chosen for all subsequent simulations, as this provided the largest physical domain size without noticeable grid effects. Additional tests were carried out to verify the phase-field simulation implementation and parameters. These tests included verifying the correct representation of the Gibbs-Thomson effect and ensuring that the interfacial energy was independent of the energy barrier coefficient. The results of the other tests to ensure the stability of the simulation are presented in Appendix A. Results and discussion 92 Figure 3.33: Grid effect in a two-dimensional phase-field simulation with two particles horizontally aligned (upper row) and slightly tilted (lower row). The tilted simulations are rotated to visualise the difference compared to the horizontally aligned simulations. In the next step, the spatial discretisation was dimensionalised to give physical meaning to the space in the simulation domain. A conversion factor α was introduced to convert the distance between two grid points into a physical length. For example, a conversion factor of 1 µm/gp would imply that for Nx = 240 gp the x-coordinate corresponds to 240 µm. The method for deriving α will be explained in the following. In this work, Nx = 240 gp and Np = 864 were selected due to feasibility and computational constraints. Under this configuration, the particles reach an average size (defined here as the sphere equivalent radius) of 14.1 gp ± 1.7 gp when occupying the space in the simulation domain with volume fractions of the particulate phase of 0.915. This value represents the average ice crystal volume fraction of maximum freeze-concentrated 10 % (w/w) sucrose and trehalose solutions, as previously determined in Chapter 3.2.2.2.2. In Chapter 1.5, published data on pore sizes in non-annealed samples were presented. These pore sizes were the product of conventional freezing and drying and were converted to circle equivalent radii of Results and discussion 93 10.0 µm to 30.0 µm [83], 12.5 µm to 25.0 µm [98], and 11.9 µm ± 5.6 µm [71]. A wide range of pore sizes can be observed within the individual publications, which can be attributed to various factors such as a heterogeneous distribution of pore sizes, differences in pore sizes depending on the axial position in the vial [98], and measurement uncertainties associated with the methods used. Unfortunately, this means that there is no uniform value for experimentally determined pore sizes that can be used as a starting point for the phase-field simulations. Therefore, a range of initial pore sizes has also been considered in this work. Due to the microscopic scale of the simulation domain, simulating a wide variety of initial particle sizes within a single simulation is impractical. Therefore, separate simulations were run for different initial particle sizes at the start of annealing. At this point, two setup options were considered:  Adapting the grid size: The first option was to derive the conversion factor α and adjust Nx (and consequently Ny and Nz) for each initial particle size. However, this approach could result in a simulation domain that is either excessively large or so small that the interfacial width would become limiting, i.e., most of the simulation domain would be occupied by interface alone if Nx becomes too small as the interfacial width is more or less fixed.  Using different conversion factors: The alternative was to vary the initial conditions using different conversion factors without altering the three-dimensional simulation domain. This meant maintaining the same number of initial particles and grid points across simulations but adjusting their physical meaning through the conversion factor. In this work, the second option was chosen in order to avoid potential excessive computational times during the execution of the simulation. The conversion factor was defined as: 𝛼 = 𝑟  𝑟  , (Eq. 20) where rexp is the experimentally determined circle equivalent radius, and rsim is the sphere equivalent radius from three-dimensional simulations. Three different initial particle sizes of 10.0 µm, 12.5 µm, and 15.0 µm were selected, based on the experimental data, and conversion factors were calculated to align the average particle Results and discussion 94 size in the three-dimensional simulation domain (14.1 gp ± 1.7 gp) with these physical sizes. For example, the conversion factor was calculated that would translate the initial physical particle size of 10.0 µm to the simulation particle size of 14.1 gp, in this case α = 0.68 µm/gp. These conversion factors were then used to calculate the base grid size for two-dimensional phase-field simulations, ensuring alignment with FDM experimental data. For a conversion factor of α = 0.68 µm/gp, the particle size at the beginning of the simulation should be 9.2 gp. This is based on the experimental finding that particles had an average radius of 6.2 µm at the earliest quantifiable time point of 10 minutes after annealing began. In order to place 100 particles in the simulation domain, each with an average radius of 9.2 gp, the base grid size was determined to be Nx = 197. In summary, this means that the conversion factors were used to determine the size of the two-dimensional simulation domains that could representatively describe 100 particles from FDM experiments. However, the various conversion factors themselves were previously calculated by comparing the established three-dimensional phase-field simulation and literature data on ice crystal sizes/pore sizes. This strategy allowed for adaptations only in the two-dimensional simulations, while the length scales of the more resource-intensive threedimensional simulations were simply adjusted by the corresponding conversion factors. A summary of all data required to calculate the conversion factors and Nx for two-dimensional simulations is provided in Table 3.5. Table 3.5: Determination of conversion factors (α) and base grid size for two-dimensional simulations (Nx2D) for the experimentally determined circle equivalent radius of particles from freeze-drying microscopy after 10 minutes of annealing at -6 °C (rexp), the initial sphere equivalent radius of particles in three-dimensional (rsim,3D), and the initial circle equivalent radius for two-dimensional simulations (rsim,2D). rexp [µm] rsim,3D [µm] rsim,2D [gp] α [µm/gp] Nx2D [gp] 6.2 10.0 9.2 0.68 197 12.5 7.4 0.85 157 15.0 6.2 1.00 134 Results and discussion 101 Table 3.6: Parameters for the calibration of the two-dimensional and three-dimensional phase-field model. The polynomial function for the initial supersaturation of the matrix phase (𝑐 ) was determined dependent on the area and volume fraction of the particulate phase (𝜑) and the base grid size of the simulation domain (Nx). 𝜑  [-] Dimensions [-] Nx [gp] Polynomial function R² [-] 𝑐   [-] 0.793 2D 197 y = 14.133x³ - 22.555x² + 12.313x – 1.3377 1 0.338 157 y = 9.9716x³ - 16.137x² + 9.0861x – 0.8157 1 0.336 134 y = 5.4224x³ - 8.8472x² + 5.226x – 0.1402 1 0.313 3D 240 y = 7.5809x3 – 13.627x2 + 8.5475x – 0.8864 1 0.369 0.828 2D 197 y = 12.225x³ - 21.028x² + 12.45x – 1.549 1 0.386 157 y = 8.1264x³ - 14.481x² + 9.0219x – 0.9637 1 0.385 134 y = 5.073x³ - 9.2468x² + 6.0656x – 0.4106 1 0.372 3D 240 y = -1.7858x3 – 0.9729x2 + 3.3107x – 0.2685 1 0.426 3.2.4.2 Coarsening parameters for the two-dimensional phase-field model In the following, the simulation parameters for the two-dimensional phase-field simulations were determined which correspond to the experimental recrystallisation rates from Chapter 3.2.2.2.3. The objective was to adjust the mobility coefficient M and relaxation coefficient L simultaneously to match the temporal evolution of particle sizes between experiments and simulations. An example of coarsening in a two-dimensional phase-field simulation is depicted in Figure 3.38. At the start of the simulation, the area fraction of the particulate phase increased rapidly until it plateaued at tsim = 300 dt, indicating that the phases had reached their respective equilibrium compositions. The height of the plateau was adjusted by the initial supersaturation of the matrix phase as discussed in Chapter 3.2.4.1. In order to reduce the total number of particles in the system, complete dissolution of individual particles was required. However, at the start of the simulation the supersaturation of the matrix phase stabilised the small particles until the equilibrium compositions were reached. Therefore, as Results and discussion 102 long as the volume flux into small particles due to supersaturation was greater than that out of them due to Ostwald ripening, the number of ice crystals did not change. After equilibrium compositions were reached, the standard deviation of the average particle area steadily increased, indicating a broadening of the particle size distribution. Consequently, some particles fell below the critical particle size and dissolved, leading to an increase in local composition within the matrix phase and slight supersaturation. In response, neighbouring particles began to grow to reduce the supersaturation. Figure 3.38: Temporal evolution of number of particles (○), area fraction of the particulate phase 𝜑 (Δ), and the average particle area (x) with standard deviation (error bars) in a twodimensional phase-field simulation with M = 1, L = 1, 𝑐  = 0.539, and Nx = 197 gp. After a lag phase of tsim = 1000 dt, the particle number started to decrease with a hyperbolic curve profile. Experimental data from freeze-drying microscopy do not show such a lag phase due to the presence of smaller particles at the time of observation (ta = 10 min), which immediately fall below the critical radius and dissolve when the sample temperature is increased. This in turn immediately reduces the number of ice crystals from the start of the observation. From tsim = 3000 dt, the decrease in the number of particles slowed down and approached an apparent plateau. The mean particle area and its standard deviation increased steadily from tsim = 1000 dt, indicating coarsening. In addition, based on the Gibbs-Thomson effect (see Appendix A), a slight increase in the area fraction of the particulate phase was observed. The simulation results are visualised in Figure 3.39. Results and discussion 103 Figure 3.39: Visualisation of the composition field in a two-dimensional phase-field simulation with M = 1, L = 1, 𝑐  = 0.539, and Nx = 197 gp. The simulation results were visualised starting from tsim = 300 dt, when the matrix phase reached its equilibrium composition. Between tsim = 300 dt and tsim = 1000 dt, the particulate phase underwent a redistribution resulting in a broader particle size distribution with the presence of both smaller and larger particles. From tsim = 1000 dt to tsim = 5000 dt, the number of particles visibly decreased, resulting in a coarsening of the particulate phase. The next step was to assign physical meaning to the temporal discretisation of the simulation. This means that the time between two iterations in the simulation (dt) had to be defined in SI units. However, dt itself is a parameter used to ensure the stability of the simulation and to minimise inaccuracies and errors. If dt is set too high, the simulation may not be able to capture the temporal evolution correctly as excessive change occurs between iterations. If dt is set too low, more iterations will have to be calculated to simulate the same phenomenon, increasing computational time. For this work, a dt of 0.05 has proven optimal, as higher values Results and discussion 104 led to crashes and lower values to prolonged simulation times. The temporal scale, i.e., what each simulation iteration means in physical time, was then adjusted based on the coarsening profile in Figure 3.38, as will be explained in the following. For M = 1 and L = 1, the estimated lag phase lasted about 1000 iterations and the number of particles remained relatively constant after 4000 iterations. In this work, M and L were chosen while the temporal conversion factor τ was determined. Technically, the temporal conversion factor could have been chosen first and then the appropriate M and L determined. Mathematically, the results would have been identical. As the experimental data covered 360 minutes of annealing and the actual coarsening occurred within simulations with M = 1 and L = 1 for about 3000 iterations, each simulation iteration was set to represent 0.1 minutes or 6 seconds. The temporal conversion factor was calculated with: 𝜏 = 6 𝑠 𝑑𝑡 , (Eq. 23) and resulted in τ = 120 s for all simulations. As previously mentioned, the recrystallisation rate in the phase-field simulations was manipulated by adjusting the values for M and L simultaneously. This approach allowed the recrystallisation rate to be modified while keeping the interfacial energy of the system constant. Altering the ratio of M and L would have affected the contributions of diffusion, dissolution and redeposition to the overall coarsening, which is outside the scope of this work. In order to investigate the relationship between M, L and the recrystallisation rate (including the temporal conversion factor), various simulations were carried out with incrementally changing M and L. The results are illustrated in Figure 3.40. It should be noted that the simulation data was not used until coarsening started to occur. This means that the starting point varied between simulations as higher values of M and L would result in a shorter lag phase. Similar to the experimental data, the simulation results displayed a linear increase in the cube of the circle equivalent diameter as coarsening commenced. The curves may appear almost stepwise. However, this is a misconception. This appearance can be explained by the fact that when particles dissolve from one iteration to the next, the total number of particles changes, leading to erratic fluctuations in the calculated average cube of circle equivalent radius. Results and discussion 105 Figure 3.40: Impact of simulation parameters M and L on the coarsening behaviour of a twodimensional phase-field simulation with 𝑐  = 0.539 and Nx = 197 gp. The slope of the regression increased with higher settings of M and L, indicating the successful increase of the recrystallisation phenomenon dependent on M and L. Following these observations, the conditions specified in Table 3.6 were used to conduct two-dimensional phase-field simulations. These simulations were aimed at determining recrystallisation rates based on the slopes of the cube of the circle equivalent radii. The results of these simulations are depicted in Figure 3.41. Figure 3.41: Impact of M and L on the recrystallisation rate in a two-dimensional phase-field simulation. Results and discussion 106 The recrystallisation rate exhibited a positive linear correlation with the simulation parameters M and L, as expected. This is because an increase in both M and L simultaneously equates mathematically to an increase in dt. This relationship can be simplified as follows In the simulations, it does not matter whether a process takes, for example, 10 minutes with a data point collected every minute, or if the process runs twice as fast with a data point collected every half minute. Ultimately, the amount of work done, and the output data remain identical. It was also observed that the recrystallisation rate scales faster with M and L as the grid size decreases. The effect of M and L on the recrystallisation rate also increased with the equilibrium concentration of the matrix phase. Using the regression lines from Figure 3.41, the simulation parameters M and L corresponding to the experimental recrystallisation rates were determined and are shown in Table 3.7. Table 3.7: Determination of simulation parameters mobility coefficient M and relaxation coefficient L dependent on experimental recrystallisation rate and base grid size (Nx). Recrystallisation rate [µm³/min] Nx [gp] Polynomial function R² [-] M & L [-] 56.80 197 y = 7.2702x + 0.5562 0.9992 7.74 157 y = 12.822x + 1.4688 0.9994 4.32 134 y = 27.29x – 1.0529 0.9999 2.12 26.13 197 y = 6.0147x + 1.5804 0.9978 4.08 157 y = 10.318x + 1.9271 0.9988 2.35 134 y = 19.78x + 1.0279 0.9988 1.27 Next, the simulation parameters M and L, as determined in this chapter, were applied to three-dimensional phase-field simulations. These simulations aimed to predict the coarsening behaviour of frozen 10 % (w/w) sucrose and trehalose solutions at annealing temperatures of -4 °C and -6 °C with initial particle radii of 10 µm, 12.5 µm, and 15 µm. Results and discussion 107 3.2.5 Advancement of the phase-field simulation to three-dimensional For the implementation of the three-dimensional phase-field simulation, all fields (e.g., composition, order parameters, etc.) were extended to include the z-coordinate. Spatial operators (∇) used for calculating gradients and divergences of these fields were accordingly modified to accommodate the additional dimension. The particle position and size distributions, as established in Chapter 3.2.3, served as the initial conditions. Volume fractions were also adjusted based on the initial supersaturations of the matrix phase, as discussed in Chapter 3.2.4.1. With the conversion factors previously determined, the following simulations are all dimensionalised in regard to space and time, i.e., volumes are expressed in μm³ and time in minutes. Figure 3.42 illustrates an example of the results from a three-dimensional phase-field simulation. Figure 3.42: Temporal evolution of number of particles (○), volume fraction of particulate phase 𝜑 (Δ), and the average particle volume (x) with standard deviation (error bars) in the three-dimensional phase-field simulation with M = 2.35, L = 2.35, 𝑐  = 0.539, and Nx = 240 gp. During the annealing simulation, the average particle size continuously increased, and it did not reach a plateau even after tsim = 360 min. Throughout the coarsening process, the standard deviation of the average particle size also increased, indicating an increasing degree of heterogeneity in particle sizes as annealing progressed. Similar to the two-dimensional simulations a lag phase was present, where the volume fraction of the particulate phase showed a rapid increase up to tsim = 10 min while the number of particles remained relatively Results and discussion 108 constant. During this initial stage, the volume fraction of the particulate phase increased due to supersaturation in the matrix phase, causing all particles to grow until reaching the desired volume fractions determined experimentally in Chapter 3.2.2.2.2. A steady increase in the volume fraction of the particles was observed between tsim = 10 minutes and tsim = 120 minutes as the simulation progressed, particularly when the number of particles decreased rapidly. As with the two-dimensional simulations, this phenomenon is based on the Gibbs-Thomson effect and is related to the curvature of the particles, i.e., the smaller a particle, the greater its curvature and the greater the variation in composition within the particle. Therefore, the volume fraction of the particle increases with coarsening. A description of this phenomenon can be found in Appendix A. The lag phase observed in the three-dimensional phase-field simulations was notably shorter compared to the two-dimensional simulations. This difference is primarily influenced by two factors. First, the addition of an extra plane of curvature in three dimensions accelerates coarsening, which effectively reduces the lag phase. Secondly, the critical particle size at which dissolution begins differs between two-dimensional and three-dimensional simulations. Specifically, the critical radius for dissolution was determined to be 6 grid points in twodimensional simulations (representing the radius of a circle) and 7 grid points in threedimensional simulations (representing the radius of a sphere), as detailed in Appendix A. Consequently, particles in three-dimensional simulations dissolve earlier and coarsen faster, resulting in a shortened lag phase. The results were also visualised and are shown in Figure 3.43. The steady coarsening of the particulate phase was clearly evident from the three-dimensional visualisations, which also showed an increase in particle size heterogeneity. Moreover, the visualisations underscored that the simulation cannot continue indefinitely. Eventually, particles could grow large enough to span from one boundary to the opposite boundary. This poses a problem as periodic boundary conditions were employed in all simulations. Under these conditions, the left border of the simulation domain effectively continues from the right border, allowing the simulation domain to be conceptually stacked indefinitely from left to right, and similarly from bottom to top and front to back, while preserving the continuity of the internal microstructure. Results and discussion 109 Figure 3.43: Visualisation of the composition field during coarsening in a three-dimensional phase-field simulation with M = 2.35, L = 2.35, 𝑐  = 0.539, and Nx = 240 gp. This also means that a particle growing from any boundary to the opposite side could potentially come into contact with itself, creating an unrealistic scenario if it becomes sufficiently large. This issue is illustrated in Figure 3.44. To address this, all subsequent evaluations were carefully monitored, and simulations were terminated if this issue arose, ensuring the integrity and realism of the simulation results. Results and discussion 110 Figure 3.44: Illustration of the periodicity of the simulation domain and the associated problems. The original simulation domain (A, red square) can be replicated and placed at each border due to its periodic boundaries. However, if the particles become too large, there is a risk that a particle will come into contact with itself (B). 3.2.6 Final adjustment of particle volume fractions In lyophilisation, the sample temperature must be reduced again after annealing to enable the matrix phase to reach its maximum freeze-concentrated state. This step prompts water molecules to precipitate into ice crystals, altering the water-solute ratio in the matrix phase, and consequently affecting its volume fraction in the lyophilisate. It is important to implement this in phase-field simulations to accurately model the particle sizes that correspond to the pores in the lyophilisate. The adjustment of the volume fractions in the three-dimensional phase-field simulations was achieved by restarting the simulation after the desired annealing step, but with an adjusted matrix equilibrium composition. Elegantly, this is also what happens physically in a real solution. As the temperature decreases, the solubility of the water in the amorphous phase changes and so does its equilibrium concentration. Figure 3.45 shows how the particle volume fraction increases as a result of supersaturation until it reaches the target volume fraction of 0.913, as calculated in Chapter 3.2.2.2.2. Results and discussion 117 Figure 3.50: Temporal evolution of the average ice crystal size in the three-dimensional phasefield simulation with different initial particle sizes (rsim,3D). The standard deviation was omitted for better clarity. It should be noted that the standard deviation became very high over the course of the simulations and has therefore been omitted from the results for clarity. This phenomenon was observed previously in Figure 3.42 and can be explained by the fact that during coarsening, larger particles tend to continue to increase in size while the rest of the particles continue to shrink until they reach the critical particle radius and dissolve. Therefore, very large and very small particles exist at the same time, which broadens the particle size distribution and increases the heterogeneity in particle size. Notably, the individual data sets do not follow a perfect curve, but tend to “jump” between annealing times. This can be explained by the complete dissolution of particles during the process. This dissolution changes the total number of particles, leading to erratic shifts in the average particle size calculations. Next, the sphere equivalent radius was calculated with: 𝑟  =  3 ∙ 𝑉  4 ∙ 𝜋  , (Eq. 25) where Vp is the dimensionalised volume of an individual particle. Results and discussion 118 It must be mentioned that the particles in this simulation were not perfectly spherical, and therefore the sphere equivalent radius is just an approximation of the particle size. However, circle or sphere equivalent radii are commonly used in the literature regarding pore sizes in lyophilisates, and facilitate the comparison with other works. The results are visualised in Figure 3.51. Figure 3.51: Temporal evolution of the sphere equivalent radius in the three-dimensional phase-field simulation with different initial particle sizes (rsim,3D) and annealing temperatures (Ta). The standard deviation was omitted for better clarity. The conversion of the ice crystal volume to the sphere equivalent radius also shows the increase over the annealing period and the dependence of the resulting particle radius on the initial particle size. The increase is not completely linear, and there may be several reasons for this. First, changes in the mass fraction in the matrix phase due to the Gibbs-Thomson effect might play a role, as larger particles exhibit lower compositional values. This phenomenon slightly influences the coarsening behaviour, as the driving force for Ostwald ripening is dependent upon the compositional difference between the particles and the matrix phase. Additionally, the interface may also contribute to the non-linear increase. Given that the interface maintains a fixed thickness, its impact is more pronounced on smaller particles than on larger ones. Consequently, small particles might appear larger in evaluations because the Results and discussion 119 interface is included in their measurement. As previously demonstrated in Figure 3.42, the standard deviation of particle size increased considerably as annealing progresses. To further illustrate this, the frequency of different particle sizes at various annealing times is depicted in Figure 3.52. Figure 3.52: Exemplary ice crystal size distributions in the three-dimensional phase-field simulation with Ta = -4 °C and rsim,3D = 10 μm. The frequency alters between the plots. Results and discussion 120 The flattening of the cumulative frequency of particle volumes, normalised between 0 and 1, clearly indicates a shift towards larger particle sizes. It is important to note that the frequencies depicted in the plots (y-axis) are not consistent across different plots. These results are qualitatively consistent across all sets of simulation parameters and confirm that coarsening in such systems leads to increased heterogeneity in particle sizes. Furthermore, it can be observed, especially in Figure 3.51, that the average particle radius at the beginning of the simulation is between 10 μm and 15 μm. This indicates that the strategy involving the conversion factor, as detailed in Chapter 3.2.3, was successfully implemented in the three-dimensional simulations and has yielded meaningful results. 3.2.7.3 Conjunction area of ice crystals The next step involved the determination the conjunction areas as defined in Chapter 3.2.7.1 and their size distributions, or more precisely area distributions. As an intermediate step, the volume resulting from the overlap of the interfaces of two adjacent ice crystals was determined, hereafter referred to as the overlapping volume. The concept of overlapping volumes is illustrated in Figure 3.53 and explained below. Figure 3.53: Concept of overlapping volume in a two-dimensional domain. The space occupied by the interface of both particles (blue and red) forms an area in two-dimensional space (green) and a volume in three-dimensional space. The overlapping volume (or area) exists because the interface of the particles exhibit a certain width, which means that some space in the simulation domain can be occupied by two (or even more) particle interfaces simultaneously. Notably, this value depends on how the interfacial width is evaluated. If the interfacial width is defined as the length between two Results and discussion 121 phases where the values deviate from the equilibrium compositions, then larger widths of 6 to 8 grid points can be measured. Alternative methods, such as in Appendix A, use the inflection point of the interface and determine lower values. Nevertheless, the thickness of the interface, given that gradient energy coefficients are not altered, remains constant. The overlapping volumes appear as thin and partially curved discs, as shown in Figure 3.54 for three-dimensional simulations. These shapes develop because the interfaces between adjacent particles deform during particle growth and coarsening processes. The total overlapping volume between all particles is given by: 𝑉 , =   𝜂  𝜂            , (Eq. 26) where Np is the total number of particles, and η is the order parameter. Figure 3.54: Three-dimensional phase-field simulation with an initial particle size of rsim,3D = 10 μm after annealing at Ta = -6 °C for ta = 360 min. Ice crystals are shown on the left, while the overlapping volumes are visualised on the right. The conjunction area that corresponds to an overlapping volume approximates about half of its surface are and can be determined with: 𝐴  = 1 2 𝐴  , (Eq. 27) where Ao is the surface area of an overlapping volume. Results and discussion 122 The concept of conjunction areas is not commonly addressed in the literature concerning the microstructure of lyophilisates. Therefore, in this work, several parameters were initially examined to identify potentially relevant aspects of conjunction areas. Similarly to how ice crystal properties were assessed, the number, average area, and total area of conjunctions were determined. An example of the temporal evolution of these properties is shown in Figure 3.55. Figure 3.55: Temporal evolution of number of conjunctions (○), average conjuncon area with standard deviation (Δ), and total conjunction area (x) during annealing in the threedimensional phase-field simulation with an initial particle size of rsim,3D = 10 μm after annealing at Ta = -4 °C. As annealing progressed, the number of conjunctions decreased. This is related to the fact that the number of ice crystals also decreased during coarsening, so that fewer ice crystals were in contact. Meanwhile, an increase in the average conjunction area was observed, which was also accompanied by an increase in its standard deviation. This indicates that, similar to the ice crystals, the system became more heterogeneous in terms of conjunction sizes. Additionally, it was noted that the total conjunction area decreased over the course of the simulation. This decline is consistent with the typical behaviour of coarsening, where surface areas within the system are minimized in favour of larger volumes. The trends in the number of conjunctions for all simulation conditions are depicted in Figure 3.56. Results and discussion 123 Figure 3.56: Comparison of various simulation parameters on the number of conjunctions in three-dimensional phase-field simulations. The number of conjunctions decreased more rapidly at higher annealing temperatures, a result of faster coarsening which reduces both the number of particles and their contacts. Furthermore, smaller initial particle sizes resulted in more conjunctions due to the increased frequency of contact among numerous smaller particles. This does not imply a change in the number of conjunctions per particle, but rather indicates that the total number of conjunctions is dependent upon the overall particle count. Next, the average conjunction area for all simulation conditions was determined with: 𝐴 󰆹  = 𝛼  ∑ 𝐴  ,       𝑁  , (Eq. 28) where α is the conversion factor, Ac is the respective conjunction area, and Na is the total number of conjunctions. This value is particularly interesting because it describes how narrow, on average, the paths for mass flow of vapour during sublimation can become. The idea is that the narrowest path acts as a bottleneck and potentially limits the total mass flow. The data illustrating this phenomenon is shown in Figure 3.57. Results and discussion 124 Figure 3.57: Comparison of the impact of various simulation parameters on the average conjunction area in three-dimensional phase-field simulations. The results demonstrate a clear relationship between the average conjunction area and the initial particle size, which is intuitive as larger particles are likely to have larger contact areas. Additionally, there appears to be a temperature-dependent influence on this relationship. This can be explained by the fact that more rapid coarsening at higher temperatures leads to the formation of larger particles and, consequently, larger conjunction areas. This temperature effect is particularly pronounced for particles with an initial radius of 15 μm, as evidenced by the steeper increases in conjunction area during annealing. Notably, the data shows erratic fluctuations, which are due to the reduction in the number of conjunctions over time. These reductions affect the calculation of the average conjunction area. Finally, the total conjunction area within the system was examined. It is important to note that the simulations accounted for different total volumes after converting to physical length units, meaning that the distance between grid points varies depending on the applied conversion factor. Therefore, it was crucial to calculate the specific conjunction area to enable meaningful comparisons of the simulation results. Results and discussion 125 For this purpose, the total volume of a simulation was calculated using: 𝑉  = 𝛼  ∙ 𝑁𝑥 ∙ 𝑁𝑦 ∙ 𝑁𝑧 , (Eq. 29) where α is the conversion factor, and Nx, Ny, and Nz are the number of grid points in x, y, and z direction, respectively. The unit conventionally used in the scientific literature to determine specific surface area is m²/g. This approach is also applicable to simulations, where the physical weight of the simulation domain can be determined by calculating its volume using Equation 29 and then applying the density of the matrix phase. The physical mass of the matrix phase in a simulation is given by: 𝑚  ,  = 𝑉  ∙  1 − 𝜑   ∙ 𝜌  , (Eq. 30) where Vsim is the physical volume of the simulation domain, ϕp is the volume fraction of the particulate phase, and ρm is the density of the matrix phase. The values for the determination of the physical mass of the matrix phase are summarised in Table 3.8. Since disaccharides such as sucrose and trehalose have similar densities, only the values for sucrose were employed in this analysis. Subsequently, the conjunction area was normalised to the specific conjunction area, defined as: 𝑆𝐶𝐴 = 𝛼  ∑ 𝐴  ,       𝑚  ,  , (Eq. 31) where α is the conversion factor, Ac is the respective conjunction area, Na is the total number of conjunctions, and mm,sim is physical mass of the matrix phase within the simulation domain. The results are depicted in Figure 3.58. Table 3.8: Conversion values for the calculation of the physical mass of the matrix phase in the simulation domain (mm,sim) from the conversion factor (α), the volume of the simulation domain (Vsim), the particle volume fraction (φp), and the matrix density (ρm). α [µm/gp] Vsim [μm³] φp [cm³/ cm³] ρm [g/cm³] mm,sim [g] 0.68 3.08 x 106 0.913 1.404[1] 3.28 x 10-7 0.85 6.01 x 106 0.913 1.404[1] 6.41 x 10-7 1.00 9.78 x 106 0.913 1.404[1] 1.04 x 10-6 1 Lescure (1995) Results and discussion 126 Figure 3.58: Comparison of the impact of various simulation parameters on the specific conjunction area in three-dimensional phase-field simulations. Generally, the specific conjunction area decreased over the course of annealing, which can be rationalised by considering the fundamental drive to minimise surface areas during coarsening. Interfaces between phases contribute to the enthalpy of the system and are consequently reduced during relaxation processes such as Ostwald ripening. Notably, the specific conjunction area decreased more during annealing when the initial ice crystals were larger and when higher annealing temperatures were used. 3.2.7.4 Surface area of the matrix phase In the next step, the surface area of the matrix phase was determined. An exemplary visualisation of the surface area is shown in Figure 3.59. As previously mentioned, it is important to note that the simulation domain corresponds to different physical volumes depending on the set of simulation parameters. This means that analogous to the specific conjunction area, the surface needed to be normalised as well to obtain comparable results between different simulations. Elegantly, the matrix surface area can be determined using values already acquired for particle surface area and conjunction area, as shown in the following. Results and discussion 133 In this thesis different initial particle sizes were used in the simulations reflecting the broad spectrum of values found in the existing literature for non-annealed samples. The results of the three-dimensional annealing simulations in the context of a frozen disaccharide solution can be summarised as follows:  Annealing leads to an increase in the average size of ice crystals and broadens their size distribution. The increase in ice crystal size during annealing, which later translates into the pore sizes in the dried lyophilisate, aligns with observations that larger pores can enhance the sublimation rate during freeze-drying.  The specific surface area of the matrix phase decreases over the course of annealing. This phenomenon might play a role in the diffusive and desorption processes during drying. Furthermore, the reduction in surface area is beneficial as it could help minimise protein agglomeration by minimising the available surfaces for protein precipitation, thereby enhancing the stability and quality of the final product.  During annealing, the number of conjunctions, which form the narrowest gaps in the porous microstructure, decreases. However, the average size of these conjunctions increases, indicating a shift towards fewer but larger channels within the microstructure. Quantifying parameters such as the conjunction area might offer better explanations why annealing effectively reduces primary drying times.  The phenomena described above depend on both the annealing temperature and the initial ice crystal size. In general, higher annealing temperatures accelerate these processes, while larger initial ice crystal sizes reduce the time required to reach plateaus. The findings from this research offer practical applications in the field of freeze-drying, specifically through the use of phase-field simulations. By simulating different annealing times and temperatures, it is possible to determine when microstructural changes reach plateau phases. This approach can help optimise freeze-drying process conditions, ensuring that the benefits of annealing (as mentioned above) are achieved efficiently. The goal is to maximise the positive impacts of annealing while minimising the duration of the annealing step, making the process more time-effective and cost-efficient. This methodology allows for a more precise control over the lyophilisation process, potentially leading to better quality products with reduced processing times. Future applications 134 4 Future applications In the following chapter some ideas for future research are explored, which can be based on the results of this work. These include both experimental and simulation topics. 4.1 Reorganisation kinetics during recalescence and crystallisation In Chapter 3.1.4, the experimental findings showed that the recalescence and crystallisation phase significantly influences the microstructure formation in a frozen solution. To reiterate, during recalescence and crystallisation, the increased mobility within the partially frozen system led to a reorganisation of the crystalline phase, resulting in the formation of ice crystals resembling grains and approximately matching the sizes of the pores within the lyophilisate. Quantifying this reorganisation event or even understanding its kinetics might be helpful for various reasons. From a scientific standpoint, it could provide an alternative explanation for the pore sizes in a lyophilisate. This work challenges the common understanding that pores in the lyophilisate directly correspond to stable nuclei formed during freezing. Instead, if the recalescence and crystallisation phase predominantly determines ice crystal size, it is plausible that multiple irregular dendritic ice crystals might coalesce to form each grain-like ice crystal observed. Thus, inferring the number of nuclei solely based on pore size observations may be misleading. From the perspective of phase-field modeling, understanding the kinetics of this reorganisation event would allow to define the starting point of the subsequent annealing simulation, since initial particle sizes would be known. However, the issue arises on how to observe and quantify such an event. If a freeze-drying microscope is used, the sample thickness of a few micrometres might already be too large for this event to occur in a pseudotwo-dimensional manner, since the initial structures seem to be much finer than the height of the sample. Resolving this challenge might require further reduction in layer thickness or potentially transitioning to a three-dimensional phase-field simulation, where the simulation height matches the distance between the two glass plates. If such a model is calibrated and validated, it would enable the prediction of pore sizes in a lyophilisate when no additional annealing step is performed. Furthermore, it would allow to optimise the freezing step in respect to the homogeneity of the lyophilisate within a vial, or when used with temperature data from the freeze-dryer the homogeneity between individual vials. Future applications 135 4.2 CFD simulations through the porous microstructure This idea could be used complementary with the simulation output from this work. The threedimensional phase-field model generates a data object that can be used as a simulation domain for CFD simulations. This integration presents a significant advantage, as fabricating such detailed and complex microstructures using alternative methods would be substantially more challenging. As mentioned in Chapter 3.2.8, the resolution of experimental methods such as μ-CT are simply insufficient to capture the microstructure , especially the thin matrix walls. The implementation of a robust CFD simulation, fed with the microstructure as a simulation domain, could be used to make predictions about the partial vapour pressure in the pores immediately at the sublimation front. The Knudsen-Langmuir equation describes the mass flow (sublimation rate) as [161]: 𝑚 󰇗 = 𝑎 ( 𝑝  − 𝑝   ) ∙      , (Eq. 34) where a is an accommodation coefficient, p0 is the saturation vapour pressure of the solid phase, 𝑝 is the partial pressure in the bulk vapour, Mv is the molecular mass of the vapour, R is the universal gas constant, and T is the interface temperature. From Equation 34 it can be derived that the purely physical sublimation rate, which does not take geometric features or product resistance into account, depends on three factors: 1) The difference between the equilibrium vapor pressure of the solid phase and the partial vapor pressure in the gas phase. 2) The temperature in the solid phase at the interface. 3) The accommodation coefficient. The Accommodation coefficient serves as a fundamental physical parameter, describing the interaction of gas or vapour molecules when they collide with the surface of a solid or liquid. The temperature at the interface, on the other hand, emerges as a multifaceted function influenced by several factors. These factors include the heat extracted by the frozen product, heat transferred from surrounding surfaces, heat released during ice crystallisation, and changes in enthalpy arising from the formation of new interfaces. Future applications 136 The pressure terms in this equation indicate that the sublimation rate is slowed down when the partial vapour pressure increases. This means that the total sublimation rate is adjusted to a balance between the pressure reduction due to migration of water vapour away from the interface and the purely physical sublimation rate that is mostly dependent on the temperature at the interface. The pressure reduction term can be assumed to be dependent on the ambient pressure, as well as pore geometry, possibly pore sizes or the smallest conjunction areas (bottleneck principle). Therefore, a CFD simulation could be used to verify if the microstructure predicted from the three-dimensional phase-field model is accurate by comparing the experimental sublimation rate with the balanced sublimation rate from the combination of CFD simulation and the porous microstructure from the phase-field model. Appendices 137 5 Appendices Appendix A – Validation and stability analysis for the phase-field simulation The following section describes multiple tests that were conducted to ensure that the simulation had been implemented correctly. On the one hand, care is taken to verify that simulation values correspond to theoretical values (see Gibbs-Thomson effect), but also that the simulation output scales properly with the input parameters (e.g., interfacial width scaling with gradient energy coefficients). i. Gibbs-Thomson (capillarity) effect on composition The composition of the phases in a phase-field model are dependent on the curvature of the interface, i.e., the more curved an interface the higher the composition of the corresponding phase. This phenomenon, namely the Gibbs-Thomson effect, is also present in the real world as the curvature-dependent increase of the saturation vapour pressure in sufficiently small particles [162]. The deviation of the composition of a particle with a curved interface is given by [163]: ∆ 𝑐  = 𝜒  ∙ 𝛾  𝑐   − 𝑐    ∙ 𝜓  , (Eq. App. 1) where χp is the mean curvature of the interface, γ is the interfacial energy density, and 𝑐  and 𝑐  are equilibrium compositions for the mand the p-phase, respectively, and ψm is a quantity defined as [163]: 𝜓  = 󰇩 𝜕  𝑓  𝜕 𝑐  󰇪    , (Eq. App. 2) where f0 is the bulk free energy density, and c is the composition. The impact of the Gibbs-Thomson effect was calculated for the two-dimensional as well as the three-dimensional model. Subsequently, simulations were performed with both models and different initial particle sizes. The results in Figure App. 1 show that simulation output and theoretical calculations agreed, and that both models exhibited a critical particle size. Once the particle radius fell below the critical particle size, the corresponding particle dissolved completely, as can be inferred from the mismatch of simulated and calculated value at 6 grid Appendices 138 points and 7 grid points for the two-dimensional and the three-dimensional simulations, respectively. Figure App. 1: Increase of particle composition dependent on particle radius for twodimensional (top) and three-dimensional (bottom) with critical particle radius (rc) for dx = 0.9. ii. Energy barrier coefficient on interfacial free energy The energy barrier coefficient Q is used in this multiphase-field model to increase the free energy at loci in the simulation domain where the interface of two adjacent particles of differing order parameters come into contact. Consequently, this penalises the overlapping of the interfaces of particles, and therefore promotes phase separation. The addition of such a coefficient allows to increase dt (increment between timesteps), because “moving” interfaces Appendices 139 due to particle growth are already inhibited by their directional expansion early on when close to other interfaces. If Q is set to 1, then more overlapping of particle interfaces is necessary to increase the bulk free energy, which in turn means that dt has to be low enough to prevent the overlapping from particles from one timestep to the next. When choosing the energy barrier coefficient, care must be taken that the interfacial free energy σ of the system is not altered, otherwise the simulation results would become dependent on Q. In order to identify suitable values for Q, a phase-field simulation was set up where interfaces were in contact with each other. Specifically, a single circular particle was placed with a thin ring of another phase separating the particle from the rest of the simulation domain, as shown in Figure App. 2, thereby creating overlapping interfaces. Figure App. 2: two-dimensional simulation domain (left) and cross-section (right) with slightly overlapping interfaces. The particle in the middle and the surrounding phase were both assigned their own order parameter to inhibit fusion. Next, the energy barrier coefficient Q was varied, and the overall interfacial free energy of the cross-section was calculated with: 𝜎 =  󰇩 𝑓 ( 𝑐 , 𝜂  ) + 𝜅   𝑑𝑐 𝑑𝑥   + 𝜅   𝑑𝜂  𝑑𝑥   𝑑 x 󰇪 , (Eq. App. 3) where f(c,ηi) is the bulk free energy density, and κc and κη are the gradient energy coefficients for composition and order parameter, respectively. The results are depicted in Figure App. 3. Appendices 140 Figure App. 3: Impact of the energy barrier coefficient (Q) on the interfacial free energy (σ) in the multiphase-field simulation used in this work. The interfacial free energy σ remained constant for values of the energy barrier coefficient Q up to 4. If Q was further increased, a small increase of σ was observed, with a subsequent steep decrease even up to negative values. The results of this test indicate that values for Q up to 4 are acceptable without altering the interfacial energy and therefore the simulation results. After some trial and error, a value of 2 was chosen for Q as it allowed to sufficiently increase dt for all two-dimensional and three-dimensional simulations. iii. Gradient energy coefficients on interfacial width and interfacial energy The gradient energy coefficients κc and κη for composition and order parameter, respectively, both influence the interfacial width ω and the interfacial energy σ in phase-field simulations. In particular, the following expressions must be fulfilled if the phase-field simulation was implemented correctly [164]: 𝜔 ~  𝜅  𝛾 , (Eq. App. 4) 𝜎 ~  𝛾 ∙ 𝜅  , (Eq. App. 5) where γ is a parameter related to the height of the energy barrier between the two phases, and κi is the respective gradient energy coefficient. Appendices 141 The interfacial energy was calculated with Equation App. 3, while the interfacial width had to be determined via the inflection point of the interface, as depicted in Figure App. 4. Figure App. 4: Determination of the interfacial width (ω) via a tangent line going through the inflection point of the composition profile. The determination of the upper and lower limits of the interface was achieved through the identification of the points of intersection between the tangent line and the extrema (maximum/minimum values) of the composition profile. Subsequently, the interfacial width was quantified by measuring the distance between these upper and lower limits. Although alternative methodologies, for example with the derivative of the composition profile, are available, this technique is commonly employed for the estimation of interfacial width due to its simple application. It is noteworthy that, in instances of interfaces exhibiting more width, this method may not encompass the entire interface. Nonetheless, the critical aspect of this procedure is to maintain consistency in the measurement. Next, and the interfacial energy and width were determined for various energy gradient coefficients, as shown in Figure App. 5. For this test, either κc or κη or both were varied simultaneously. The results show that the proportionality expressions from Equation App. 4 and Equation App. 5 were fulfilled, as the interfacial width as well as the interfacial energy correlated linearly with the square root of the energy gradient coefficients. Appendices 142 Figure App. 5: Dependency of interfacial width (ω) and interfacial energy (σ) on the gradient energy coefficients with γ = 1. The interfacial width depends on the evaluation method.