Numerical Modelling of Formation of Highly Ordered Structured Micro- and Nanoparticles – A Review †

The aerosol particles play a significant role in the environment and human health. They are also increasingly used in medicine (drug carriers), preparation (nanocatalysts) and many other fields. For these applications, the particles have to possess unique properties which arise directly from their structure and topology. Indeed, the functionality of the nanostructure particle is defined through its application, like chromatography, sensors, microelectronics, catalysis, and others. That is the reason why people are more and more interested in manufacturing structured particles. The structured particles are the particles with well-defined topological structure. Examples of such particles are porous particles, hollow particles (with the empty space inside), or multi-component particles with the segregation of components in the particle structure. Such particles usually have very interesting features, e.g. porous particles have a significantly larger surface area than the simple spherical particles with similar volume. The present paper contains a comprehensive review of the numerical simulation methods of the formation of highly ordered structured particles. The most important methods will be described in detail and their fields of application (with specific examples), advantages, limitations and information about their accuracy will be given.


Introduction
The structured particles are the particles with welldefined, usually ordered, topological structure (Yang S.-M. et al., 2008). Examples of such particles are porous particles, hollow particles (with the empty space inside) or the multi-component particles with the segregation of components in the particle structure (e.g. core-shell particles).
Such particles usually have very interesting features, e.g. porous particles have a significantly bigger surface than the simple spherical particles with similar volume. This feature allows us to use these particles as a nanocatalytic objects. Even if the substance that builds the particle has not catalytic properties, the particles may be covered with the thin layer of catalyzing substance.
Hollow particles (HP), on the other side, may be considered as a carrier for encapsulation of molecules or nanoparticles (Liu B. and Hu X., 2020). Other possible applications of HP are coatings (polymeric HP are considered as a white pigment, competitive to conventionally applied TiO 2 ), cosmetics (carboxymethyl cellulose-based HP are applied for encapsulation of active ingredients of various cosmetics) and electronics (mesoporous polymeric HP are considered as a filler in the electrolyte membranes for solid-state batteries). For the possible application of HP see the comprehensive review of Wichaita et al. (2019) and references therein.
Core-shell particles (CSPs) are a class of particles which contain a core region and a shell layer covering a core. The core and the shell can be different materials or the same materials with different structures (e.g. solid core and highly porous shell). One of the possible applications of core-shell particles is chromatography, e.g. high-performance liquid chromatography, where silica core-shell microspheres are used as a packing material (Hayes et al., 2014).
CSPs are an example of composite particles. In general, by the concept of composite particles we understand particles made of at least two different substances. The atoms or molecules that make up a particle may be randomly distributed in its structure or create some geometrically ordered pattern. In the latter case, it can be said that the composite particle is a structural particle. Apart from CSPs, an example of these particles are the so-called Janus particles (JPs). They were named after the Roman god Janus, who had two faces, each facing a different direction. Similarly, Janus particles have two sides which exhibit different properties. E.g. one of these sides may be non-polar and the other polar. Such particles may exhibit properties and behavior similar to surfactant molecules (de Gennes, 1992). Another example of Janus particles is when both sides are made from different materials (e.g. different metals), which may catalyze different chemical reactions. Such particles may catalyze a chemical reaction and, at the same time, move to regions where the concentration of reactants is higher.
One should also not forget about nano-and microcrystals. Such small crystals are often monocrystals i.e. the crystal lattice of the entire particle is continuous and has no defects nor grain boundaries (however, there are exceptions). Thus, the structure and morphology of nano-and microcrystals are defined by the geometry of crystal lattice. One of the most known examples are the snowflakes. Though they may differ in appearance and shape ("no two snowflakes are alike"-as the famous saying goes) but all of them have characteristic six-fold symmetry which arises from the hexagonal crystalline structure of ice.
The structure of matter, on a microscopic scale may dramatically influence its functionality. Engineering of particle structures on a microscopic scale would also lead to their unique optical, electrical, magnetic, or rheological properties maintaining at the same time the high surface area. The functionality of the nanostructure particle is defined through its application, like chromatography, sensors, microelectronics, catalysis, and others.
Well-controlled synthesis of nanomaterials and nanoscale characterization enables us to unambiguously correlate the structural properties with the physical, chemical and biological properties of nanomaterials, which form the core of nanoscience research. Excellent examples of the use of advanced particle engineering to design the structure and properties of catalysts are described in the review paper (Czelej et al., 2021). An essential requirement for nanomaterial synthesis is the uniformity of the final product in size, shape and chemical composition. Recently, many synthesis approaches have been developed to produce high-quality nanoparticles, nanorods, nanowires, or other nanostructures using metals, semiconductors and oxides (Mao et al., 2005;Mohanty, 2011).
Each of the types of structural particles mentioned can usually be obtained in several different ways. These methods depend on the particle structure itself and will be briefly presented in the next chapter. There is a suitable synthesis method to obtain a particle with a predetermined structure. A good guideline that helps to choose the appropriate method of particle synthesis may be to perform a numerical synthesis of such a particle. Numerical calculations, with the participation of a computer, allow for cheap and quick prediction of the effect of the assumed procedure. For this reason, in parallel with the methods of synthesizing particles of a given structure, methods of numerical modeling of the synthesis of such particles are being developed. Therefore, the chemical development of numerical methods related to the preparation of structural particles has taken place in the last 30 years. Currently, there are many methods that allow obtaining particles with a specific structure.
The purpose of this publication is to provide an overview of the numerical methods used to describe the formation of particles with a specific structure. Four of them we can call classic: these are groups of methods based on heat and mass balances, population balances, single atoms or molecules tracking (based on molecular dynamics) and Monte Carlo methods. They have been applied many times in numerous practical issues and constantly developed. Apart from them there will be described some less popular methods which may still be useful in specific issues.

Types of structured particles and method of their synthesis
Currently, many methods are known for the synthesis of structured particles. To choose the right one for a specific process, you must first determine what type of particles you want to get. Methods that will allow synthesizing, for example, HP, may not work at all in the case of JP and vice versa.
In general, all the methods of formation of particles, independently -structured or not -may be divided into two main groups: top-down and bottom-up (Verma et al., 2009;Fernandez et al., 2014). The first group of synthesis methods implies that the nano-or microparticles (NMPs) are synthesized by breaking larger structures (e.g. crushing, grinding or etching out crystals). A top-down method can thus be considered as an approach where the building blocks are removed from the substrate to form the particles (Aryal et al., 2019).
Oppositely, a bottom-up synthesis method means that the particles are synthesized by joining smaller structures together (e.g. atoms, molecules or smaller nanoparticles). A bottom-up approach can thus be viewed as a synthesis approach where the building blocks are added onto the substrate to form the particles (Betke and Kickelbick, 2014;Lesov et al., 2018).
From a practical point of view, the bottom-up synthesis is considered more advantageous than the top-down one, because the former usually leads to more homogeneous products with less defects and more homogenous chemical composition. Thus, in the next part of this paper, primarily the bottom-up methods will be considered.
Ones of the earliest known methods of synthesis of NMPs were the crystallization of particles from solutions (and similar precipitation) and spray drying (de Souza Lima et al., 2020). These methods have led usually to poorly structured particles. Indeed, by means of above-mentioned techniques particles with a wide range of sizes and shapes and a chemical composition that is difficult to control are obtained. Nevertheless, it is possible to obtain structured particles by means of these methods under strict conditions. Crystallization and precipitation may lead to the formation of nanocrystals and spray drying-to spherical particles, or even core-shell or hollow particles.
However, the more complex particle morphology, the more advanced methods should be used. Template methods may serve as an example of such a sophisticated and multistep approach, where first the template is synthesized, then modified to obtain the particle, and finally removed to create desired pores and voids in the structure. E.g. in the case of HP formation, the template is a simple spherical particle. This particle is then coated with a layer of the desired material and finally -removed, e.g. by means of etching (Liu B. and Hu X., 2020), leaving an empty shell. Following this description, we may conclude that the template method may be considered a combination of bottom-up (template formation, coating) and top-down (template removal) methods.

Hollow particles
Hollow particles (HPs), defined as particles (usually of spherical shape) with void space surrounded by a shell are ones of the most popular structured particles (Lou et al., 2008;Wang X. et al., 2016;Liu B. and Hu X., 2020). It has been observed that this type of particle morphology may appear during the spray drying technique when the solution or slurry is atomized into small droplets which are then evaporated (Mezhericher et al., 2010;Vehring, 2008). This technique, together with a similar gas blowing technique, was the most popular synthesis method of this kind of HPs (Kellerman, 2011;Tsotsas, 2012;Liu B. and Hu X., 2020).
Today, a lot of other methods for the formation of HPs are known. The most popular ones are collectively referred to as "colloidal templating" (Caruso et al., 1998). This group of methods boils down to producing a template on the basis of which the shell of HP is synthesized, and then removing this template (Tu and Gösele, 2005;Fan et al., 2007). In the group of templating methods we may distinguish hard-templating methods, soft templating methods and self-template methods (Liu B. and Hu X., 2020). During the hard templating methods hard templates with specific (usually spherical) shapes are synthesized first, and then their outer surface is coated with a layer of the desired material. During the soft-templating method, the role of "template" plays the droplet of emulsion or the bubble of gas (Hadiko et al., 2005). Finally, the core material is selectively removed to obtain the hollow structure (Lou et al., 2008). Finally, self-template methods denote direct synthesis without the need of additional templates (Zhang Q. et al., 2009). Methods from this group contain two stages: the synthesis of template nanomaterials and the transformation of templates into HP. The main difference between the conventional templating methods and self-templating methods, the templates used in self-templating methods are not only the templates used to formation of a void space inside the HP but also form the outer region -the crust of HP.
An exhaustive review of methods of synthesis HPs is presented in (Liu B. and Hu X., 2020).

Core-shell particles
Core-shell particles are the particles which contain an inner core and an outer shell (Ghosh Chaudhuri and Paria, 2012). The core and the shell differ in material type, in the structure (e.g. solid core and porous shell made from the same substance) or both of them (Hayes et al., 2014).
Core-shell particles are usually synthesized by a multiple-step process. The core is synthesized first and the shell is then formed on the core particle via different methods, depending on the type of core and shell materials and their morphologies (Ghosh Chaudhuri and Paria, 2012). Among these methods, layer-by-layer approach via electrostatic interaction, shell synthesis on pre-formed cores, droplet-based microfluidic method and many others can be mentioned. An exhaustive review of methods of synthesis of core-shell particles has been made in work (Hayes et al., 2014).

Janus particles and other composite particles
As it has been stated in Introduction, Janus particles (JP) are the particles which have two sides which exhibit different properties. They are probably the best-known and most widely used composite NMPs. The excellent overview of these particles' properties, applications and synthesis methods is reviewed in the paper (Walther and Müller, 2013).
There are three principal methods of synthesizing JP (Walther and Müller, 2013). The first one is based on the self-assembly of block copolymers (Walther et al., 2007). As a result of this process, individual blocks of polymers form regions of a particle with different properties. This method allows obtaining JPs with various shapes, not only spherical ones (Liu Y. et al., 2003). The second method uses "masking" (blocking) of a part of the surface of an already formed particle and then surface modification of the unblocked part of the surface (Lin et al., 2010). After removing the masking agent, a particle is obtained having surface regions with different properties. Finally, the third method consists in creating droplets of emulsions composed of various substances, and then leading to phase separation in the droplets thus formed and their solidification. In addition to JPs, there are many particles that deserve to be called composite particles because of their complex, usually more or less ordered internal structure containing at least two substances with different properties. Some of them can be treated as certain modifications of the more popular structures described above. For example, if in a core-shell particle the outer shell is not a continuous, solid or porous spherical layer but a layer of smaller particles covering the core more or less closely, we are dealing with a raspberry-like particle (when the core is tightly covered with smaller particles) or core-satellite particle (when there are larger gaps between the particles covering the core and possibly the particles covering the core are of different sizes) (Mourdikoudis et al., 2015). Particles analogous to HPs but with multilayered crust are sometimes reported as "sandwich-like particles" (Liu Q. et al., 2020). Yolkshell particles, like CSPs, have an inner spherical core and an outer spherical shell -but in their case, the shell is "too big" for the core: a small spherical nucleus can move freely in the space surrounded by the shell (Sun X. et al., 2020). There are many more such examples and it is to be expected that as nanotechnology develops, structured particles with an increasingly complex internal structure will be obtained.
These particles are usually synthetized by means of gasphase or vapor-phase methods (Swihart, 2003) or by means of coating of simple particles with appropriate substance (Demirörs et al., 2009). Very promising method of synthesis of composite particles is flame spray pyrolysis (FSP) (Mädler et al., 2002;Chiarello et al., 2008). In FSP, the solution containing liquid precursor and fuel is atomized to form the droplets which are then evaporated and ignited by means of an additive small flame (so-called pilot flame) to form solid nano particles.

Other types of structured particles
In addition to the particles with the structure described above, other, less popular types can also be found in nanoengineering applications. Among them, porous particles, as well as superficially porous particles, are worth mentioning. The latter kind of particles has a structure somewhat similar to CSPs, where the core is a solid sphere, while the shell is characterized by the presence of pores of various sizes. Both of these particle segments can be made of the same substance and differ only in terms of internal structure-so it is not a composite particle. Porous and superficially porous particles are applied as catalyst supports (Ye et al., 2017) and filling in chromatographic columns (Wang X. et al., 2012).
Micro-and nanocrystals of various substances also may be considered structured particles. They have a welldefined shape which in turn results from the structure and symmetry of the crystal lattice of the substance under consideration. The formation of submicron-sized single crystals has long been studied both experimentally and theoretically.
Last but not least, viruses should be mentioned as an example of structured particles. These objects, which are intermediate between particles and living organisms, are characterized by a ordered structure, which consists of the nucleic acid (DNA or RNA), a shell (capsid) made of proteins and sometimes some other additional components. The structure of viruses is much more complicated than the structural particles discussed earlier, e.g. CSPs. Due to the great importance of viruses in medicine and biomedical engineering, there are recently more and more attempts to describe the formation of this type of "particles" numerically (Kraft et al., 2012). Most attention is paid to the modeling of capsid formation from protein subunits (Zlotnick, 1994).
The most popular and most often investigated types of structured particles are presented in Fig. 1.

Models based on heat and mass balance equations
One of the most popular groups of numerical models of the formation of NMPs are those based directly on equations of heat and mass balances (HMB). They are especially popular for considering droplet formation during the spray drying method; however, they are also applied to model some other processes of formation of NMPs, e.g. precipitation. Fig. 1 Examples of structured particles: (a) hollow particle, (b) coreshell particle, (c) Janus particle, (d) raspberry-like (core-satellite) particle, (e) yolk-shell particle, (f) superficially porous particle.
The aim of all those models is to formulate the mass balance equation which takes into account the mass flux from surrounding to the particle or vice versa. In parallel, the energy balance equation is formulated. Both of those equations, with proper initial and boundary conditions are solved using widely applied numerical schemes, e.g. finite difference method. Levi-Hevroni et al. (1995) and Mezhericher et al. (2007;2008; have applied this scheme for the spray drying process. During the first stage of drying-evaporation of droplet-they considered the escape of solvent molecules from the droplet leading to a decrease in its radius (Fig. 2a). The mass balance equation -when assuming spherical symmetry of a droplet-takes a form of the ordinary differential equation describing the dependence of droplet radius on time (Levi-Hevroni et al., 1995): Where the mass flux ṁ v is given as (Mezhericher M. et al., 2007): All the symbols are described in Nomenclature at the end of the paper.
The mass transfer coefficient h D is obtained usually from the Ranz-Marshall correlation (Levi-Hevroni D. et al., 1995) and is a function of the diffusion coefficient of solvent vapor in air as well as viscosity and density of air. These parameters, together with solvent density ρ d and saturated solvent vapor concentration ρ v,s , are dependent on the temperature of a droplet surface which may be obtained by solving a heat balance equation: With the boundary condition at the droplet surface: The set of equation (1-4) is solved through the whole first drying period by means of well-known discretization methods (e.g. finite difference scheme). This period ends when the concentration of solutes or suspended particles in a droplet is big enough to start the formation of solid crust. The simulation of the second drying stage is much more complicated. In their early works, Mezhericher M. et al. (2007;2008) considered the emerging crust as a porous solid medium through which further solvent evaporation took place. Because the further decrease of particle radius is impossible after the formation of the crust, the evaporation of a solvent results in formation of internal bubble inside the dried particle ( Fig. 2b). That leads directly to the HP formation.
Later on, these authors have considered the more realistic crust structure as built from spherical nanoparticles suspended in droplet (Mezhericher et al., 2011;2012). This approach allowed formulating the criterion of the final particle morphology: hollow particle or solid particle. Similar investigations -with the use of the same approach -have been done to investigate the morphology of particles obtained by means of spray pyrolysis (Lenggoro et al., 2000;Eslamian et al., 2006). The described model allowed the formulation of the conditions of solid or hollow particle formation for this process as well. Also, Shabde et al. (2005) have used the described method to simulate the formation of HP during spray drying. Likewise, Maurice et al. (2017) considered the dynamics of growth of internal bubble during the HP formation in spray drying process.
The main difficulty in the application of the numerical method described above is the existence of moving boundaries. Indeed, during the drying a droplet shrinks and thus the condition (4) is fulfilled for the different value of r at every time step. The Mezhericher group usually applies the method elaborated by Moyano and Scarpettini (2000) for one-dimensional moving boundary problem with implicit boundary conditions. They use the Landau's transformation in order to work with a fixed number of nodes at each time step despite the moving boundaries. Shabde et al. (2005), on the other hand, have used the gradient-weighted moving finite element method (Miller, 1997). Recently, a work (Larbi et al., 2021) appeared in which a spectral collocation method has been applied. This method is based on the expansion of the unknown function in terms of the Lagrange elementary polynomials. The solution of HMB equations comes down to finding the dependence of the coefficients of this expansion on time. Another step towards the description of the structure of particle formed in the spray drying process is analyzing the distribution of concentration of solutes inside the droplet. The mass balance equation for the solute takes the form of convection and diffusion equation: The radial velocity of a solvent in a droplet v r is dependent on the velocity of the droplet surface caused by its shrinking. Eqn. (5) may be transformed into: where r d is a radius of a droplet and d R r r  denotes the normalized radial coordinate. Eqn. (6) is a base for the numerical model of particle formation elaborated by Vehring et al. (2007) and further developed by other authors under the name of VLF model (Chen X.D., 2011). Originally, the model assumed steady-state evaporation of a droplet what defined a proper boundary condition at the droplet surface. Later on, it has been extended to take into account the transient nature of the solute concentration profiles (Boraey and Vehring, 2014).
As a method of description of particle formation during spray drying, the VLF model has been applied mainly for the simulation of HPs formation. It allows finding out the time of drying of a droplet and the final thickness and density of a crust.
More advanced approaches use computational fluid dynamics (CFD) methods to solve HMB equations, e.g. finite element method (Vasiliauskas et al., 2015). Recently, Wu and Chung (2020) used the volume of fluid method to simulate the synthesis of HPs by means of spray drying process. At the expense of a significant complication of the mathematical apparatus, this approach makes it possible to follow the dynamics of particle synthesis without introducing additional assumptions, e.g. regarding the symmetry of the particle. The proposed model allowed directly distinguishing four sequential stages of the drying process: temperature adjustment, isothermal drying, crust formation, and crust thickening. While in the previously discussed papers, these stages were artificially separated, and for each of them, a slightly modified system of equations was used-in this case, one system of equations allows describing the entire drying process.
The methods based on mass and heat balance equations have also been applied for modeling of synthesis of NMPs other than HPs, e.g. core-shell particles. The works presenting the results of this type of simulation focus not on balance within a single emerging particle but rather on the entire system (apparatus) in which the particles are formed. In such an approach, the material balance, which describes rather chemical reactions taking place during the process than the convection and diffusion of compounds, is coupled with-usually simplified-population balance of particles, e.g. (Aul'chenko S. and Kartaev E., 2018). The authors use a single-liquid flow model to describe the synthesis of composite core-shell TiO 2 -SiO 2 particles from the gas substrates (mixture of titanium tetrachloride and silica tetrachloride). The efficiency of solid particle formation is obtained from kinetic reactions-both homogeneous and heterogeneous-and phase transition between gas and solid phase. Such an approach, coupled with a simplified balance of particle number, allowed to calculate the average size of particles at an arbitrary time point. From this model any other properties of structured particles could be obtained, e.g. the mean thickness of silica shell-however, authors did not consider this properly.
Model based on heat and mass balances equations has been also applied for simulation of FSP method (Ren et al., 2021). Similarly as in work by Mezhericher group, authors of this works considered single droplet and assumed its spherical symmetry. The main difference between modeling of spray drying and FSP is considering the effects of chemical reactions in the last one.
A similar assumption has been applied in (Hamzehlou et al., 2016) to simulate the dynamic evolution of the morphology of polymer-polymer latex particles. This paper's main goal was to consider the formation of clusters of one of polymers (P1) at the surface of the particles built from the second polymer (P2). By considering the distribution of cluster size at different time points, the equilibrium morphology of particles (described as core-shell or hemispherical particles) and non-equilibrium one (a few of P1-clusters with various sizes deposited randomly on particle P2surface) were recognized.
To sum up-a group of methods based on mass and energy balance equations has been found and still is widely used to describe the synthesis of particles with various morphologies. It also allows tracing the dynamics of the particle morphology over time and obtaining information about the quantitative distribution of parameters describing this morphology. With the help of this method, it is also possible to determine the structure of non-equilibrium particles (which are sometimes formed under not properly selected process conditions). The most important disadvantage of the described method is its relative complexity and high demand for computing power-and for this reason, other numerical methods are constructed, which sometimes constitute a better choice when it comes to particle synthesis description. These methods are described in the next section.

Multiple shells model
As mentioned above, the most challenging issue concerning the solution of differential balance equations (e.g. modeling the droplet during spray drying process) are the moving boundaries of drying droplets. To get around this difficulty, a set of methods in which the droplet is divided into imaginary regions were proposed. These regions are characterized by the same value of the concentration of every solute and the same temperature. This is, of course, an idealization of real situation and these regions mimic the parts of dried droplet where these parameters are nearly the same. This group of numerical methods is usually called the distributed parameter methods (Wang S. and Langrish T.A.G., 2009). Because of the spherical symmetry of the drying droplet, the above-mentioned regions usually have a shape of concentric spherical shells (the central one is a full sphere); thus, the described method is also known under the name multiple shells model (Parienta D. et al., 2011).
The first method from this group simply utilized the mass balance equations discretized in respect to radial coordinate (Adhikari B. et al., 2007). Wang S. and Langrish T.A.G. (2009) have used an explicit form of the concentric shells to analyze the changes of solute concentration and temperature in every part of a droplet. Parienta D. et al. (2011) used the same model for investigations of diameter and composition change of exhaled respiratory droplets. Also, similar model has been applied to investigate the distribution of protein and sugar in an evaporating droplet (Grasmeijer N. et al., 2016a;2016b). In the last work not only the droplet but also the layer of surrounding air was divided into concentric shells. This allowed the solvent evaporation rate from the droplet to be simulated directly rather than using mass transfer coefficients obtained from semi-empirical correlations.
All these works, however, did not describe the process of solid particle formation. In fact, as they analyzed the change of solute concentrations inside the droplets, their results may be accounted as numerical modeling of the first step of particle formation.
The second step, i.e. crust formation, has been modeled in (Gac and Gradoń, 2013a). The scheme of the method has been presented in Fig. 3. The authors assumed that the crystallization in this shell starts and the solid crust forms when in any of the shells (typically, it was the outermost shell) the supersaturation exceeds the critical value. In forthcoming time steps, this shell is not considered -the second shell becomes the outermost. During this period, the diameter of all the shells stops decreasing and remains constant till the end of the simulation. Instead, the concentration of all the compounds (including a solvent) in the central sphere decrease and finally reaches zero. When it is happened, in this sphere there is no solvent nor solutes -that means the gas bubble is formed. In next steps the concentration of solutes in the outermost shell increases (what may finally lead to the attachment of this shell to the solid crust) and the concentration of both solutes and solvent in the innermost shell decreases (what may result in the attachment of this shell to the central bubble). This scheme is repeated until all the shells do not belong to solid crust or gas bubble.
The multiple shells model, as described above, has found its application not only for modeling the formation of NMPs. There are a few works which use a similar model to describe the thermal evolution of planetesimals (Weiss and Elkins-Tanton, 2013;Blackburn et al., 2017). These works describe the changes of temperature profile inside the spherical planetesimal and the change of aggregation state of its compounds -they, however, do not consider the change of the diameter of the body. In (Gac and Gradoń, 2013b) the multiple shells model has been adapted to simulate the evolution of deposits appearing during the drying of sessile droplets (e.g. the famous "coffee rings").
It should be stated that the multiple shells model is characterized by a very simple mathematical apparatus. The transport of mass and energy is described here entirely using ordinary differential equations, which can be solved with not very advanced numerical techniques. However, for such a description to be applied, the simulation object (the resulting particle) must have an appropriate (usually spherical) symmetry. This condition, in principle, limits the scope of potential applications of the described model to hollow particles or core-shell particles. However, it is such a common morphology of particles produced in various processes that the described model can be considered, despite significant limitations in its applicability.

Models based on population balance equations
Population balance equations (PBEs) are widely used in several branches of modern science, mainly in chemical engineering, to describe the evolution of a population of particles. They define how populations of particles develop specific properties over time, taking into account such phenomena as nucleation, aggregation, coagulation, breakup of aggregates, evaporation, dissolution, etc. (Vetter et al.,Fig. 3 Scheme of the multiple shell model of (a) a dried droplet and (b) a wet particle with solid crust (marked in black and white), wet core and a gas bubble inside. 2013; Hao et al., 2013).
Moving on to the modeling of synthesis of NMPs, the PBEs may be applied to describe the evolution of nanoparticles' agglomeration leading to the creation of the larger structured particle. This procedure is usually applied to model the formation of particles by means of spray drying. Handscomb et al. (2009a) applied the PBE method for modeling spray drying kinetics. As during spray drying HPs are often formed, the authors investigated the time it takes to form a solid crust for various conditions of the process. For a description of the drying of a droplet containing suspended solids, the dynamics of particle number density function inside the droplet N(L,r,t) is analyzed, where L is particle size, r-its radial coordinate (spherical symmetry of dried droplet is assumed) and t-time. The equation fulfilled by function N has a form: where G is the velocity of motion through the state space of internal coordinates, i.e. the linear growth rate. The authors have assumed that this growth rate is independent of particle size. For solving Eqn. (7) a moment method is employed to yield the evolution of the first four integer moments of the internal coordinate. In their subsequent work (Handscomb C.S. et al., 2009b), the same authors extended their model for droplets drying in the presence of a solid crust. The improved model incorporates crust thickening, along with wet and dry crust drying. That allows predicting not only the kinetics of a particle formation process but also its final structure: crust thickness and its porosity. The PBEs method of modeling the formation of particles during the spray drying process has also been applied in (Bück et al., 2012). The authors considered the evolution of the distribution of nanoparticles suspended in dried droplets with particular emphasis on their aggregation. In contrast to the above-mentioned works, the effect of droplet shrinking on the PBE was not modeled as an artificial boundary source but as a flux term in the PBE. The authors of the discussed study also drew attention to the hierarchical structure of the solid crust formed during drying. They defined its porosity as ε sh , taking into account that the crust was formed from previously formed aggregates which also have internal porosity ε agg . The model proposed by the authors allowed finding out both the porosities. However, the potential of these results has not yet been fully exploited.
The modeling of synthesis of structured particles by means of PBEs seems to be one of the more general methods. A significant feature of this approach is a lack of detailed assumptions about the morphology of the formed particle. E.g. one need not assume that the final particle is of a spherical shape or spherical symmetry -these proper-ties are to be obtained directly from the model. The downside of the described method may be a relatively large complication of the mathematical apparatus. However, modern computers make it possible to overcome this difficulty.

Models based on single atoms, molecules or nanoparticles tracking
As it has been mentioned in Sec. 2, many structured NMPs are synthesized by means of bottom-up methods, i.e. as a result of joining of atoms, molecules or smaller particles (hereinafter referred to as the common term "individuals") into the desired structure. Therefore, an obvious method of numerical modeling of such a process is a direct tracking of the mentioned individuals making up the particle. This may be done by numerical solution of the equations of motion of individuals in form: In Eqn. (8) ẍ denotes the second derivative of a position vector of individuals, i.e. the vector of acceleration and F is a vector of force, which may (or may not) depend on position vector, velocity vector and time. The most common method in this group is molecular dynamics (MD) which is widely used for the description of various phenomena occurring at the atomic or molecular level.
MD algorithms are common for the simulation of the behavior of structured particles at various conditions. For example, the alloying reaction in core-shell Al-Ni particle (Levchenko et al., 2010), changes in structure during the heating and cooling of such a particle (Song and Wen, 2010), melting behavior of core-shell Pt−Au nanoparticles ) and many others have been investigated by means of this method.
Of course, MD has also been applied for modeling of synthesis of structured NMPs. Förster et al. (2019) used a very simple MD model to investigate the dependence of a particle structure on the parameters of the interactions between the atoms. In their investigations, three kinds of atoms have been considered. The first kind represented an inert gas (the atoms of this species react neither with each other nor with atoms of other kinds). Two other species are reactive and interact with each other by means of classical Lennard-Jones (L-J) potential: In Eqn. (9) i and j denote the name of species of atoms (referred as A and B), ε ij is a parameter describing the strength of interaction and σ ij -a range of this interaction. By changes of these two parameters, three types of structures of NMPs have been observed. The most important parameter influencing the morphology of the obtained particles turned out to be ε ij . It appeared that when ε BB < ε AA < ε AB (i.e. the interaction between atom A and atom B is stronger than both between two atoms A and two atoms B) an alloy particle is obtained. When ε BB < ε AB < ε AA , most of the simulations ultimately result in a core-shell particle. When ε AB < ε BB < ε AA , the Janus particle is most frequently obtained, however also quite often (about a third of the cases), the formation of pure B particles and pure A particles is observed.
Another example of the application of MD for the formation of structured particles has been applied in (Singh et al., 2014). In this paper, an MD was used to model the gas-phase synthesis and growth of silicon-silver (Si-Ag) hybrid nanoparticles. For simulating such a real system more advanced model of interatomic interactions has been applied than L-J potential. While the latter takes into account the interactions of pairs of atoms, it may be applicable with good accuracy to describing rarefied gases, especially noble ones. However, this model may be too vague to describe, for example, metals in solid state. In this case, L-J potential allows insight into qualitative properties of an inner structure of formed particles. More specific problems require taking into account collective effects between atoms. One of the most popular models including these effects is the embedded atom method (Foiles et al., 1986;Daw et al., 1993) which is a semi-empirical, many-atom approximation describing potential for computing the total energy of a metallic system. This model has been applied by Singh et al. (2014) to describe the energy of interactions between Si and Ag atoms. The nanoparticles obtained in this process had a form of Janus or core-satellite particles. The simulations allowed concluding that the individual species first form independent clusters of Si and Ag without significant intermixing. The formation of stable hybrid structures took place in the subsequent stages of process. In addition to MD, other models are used in engineering applications, which also boil down to tracking the evolution of the system by solving equations describing the dynamics of individual particles included in it. An example of such model is dissipative particle dynamics (DPD)-a coarsegrained model, in which functional groups or polymer segments are considered as individual "particles" (coarsegrained beads) which interact with each other (Español and Warren, 1995). The forces acting between beds i and j are given by the formula: where r ij is a distance between beds i and j, r c is a range of interactions between beds and ε ij -repulsion parameter, characterizing the magnitude of the force of interaction.
The DPD method has been applied in (Chen H. and Ruckenstein, 2012) to investigate the formation of complex colloidal particles formed by self-assembling a three-arm star-shaped polymer. By using this mathematical apparatus, pathways for a few kinds of structured particle formation have been identified, e.g. core-shell, multicore or stacked layered particle. The control parameters in the model were repulsion parameters between polymer segments and between polymer segments and solvent particles.
Similar methods can be used to track the movement of nanoparticles as they are used to create larger particles. Jabłczyńska et al. (2018) presented an experimental and modeling study of the formation of structured microparticles by spray drying of a droplet of slurry containing nanoparticles in one or two sizes. The change of droplet radius during the first stage of drying has been described by means of mass and heat balance equations as described in Sec. 3.1. The motion of nanoparticles inside the droplets has been tracked by direct solving of their equations of motion. To take into account the diffusion of nanoparticles inside the dried droplet, these equations have a form of overdamped Langevin equations: Here b denotes a damping coefficient and Γ(t) is a stochastic component of the motion. This component represents a white noise with zeroth mean value, i.e. it fulfills the conditions: It has been shown that as a result of spray drying of such droplets, solid (full) aggregates or hollow particles built from nanoparticles may appear. Implemented scenario depends on the process conditions: the temperature of the drying gas and the initial diameter of the slurry droplets. In the case of two sizes of primary nanoparticles and the formation of HP, self-organization of nanopraticles inside the crust may also appear. The larger nanoparticles form a hexagonal pattern while the smaller ones are trapped in void spaces between larger ones.
The molecular dynamics and related methods should be considered by far the most accurate models of the formation of structural particles. Indeed, in the case of these techniques, we directly follow the behavior of atoms, molecules or nanoparticles-i.e. the "building blocks" that form an organized particle. The main disadvantage of this group of methods is the long computation time. Especially in the case of the original MD method, a very large number of differential equations should be numerically solved, corresponding to the number of molecules or atoms involved in the particle synthesis. It should be emphasized that even a small particle consists of many hundreds of thousands or even millions of atoms or molecules (e.g. Förster et al., (2019) conducted simulations for 200,000 atoms). Another problem may be to formulate a realistic model of interactions between atoms and molecules that make up the particle. Especially when dealing with a composite particle composed of metal atoms, one should take into account the complicated effects resulting from the quantum mechanics of many bodies. However, it is worth paying attention to the fact that even a greatly simplified description using the L-J potential allows obtaining qualitatively (though not necessarily quantitatively) correct results.
The methods from the group discussed in this section should be used when interested in a detailed description of the particle structure. The models presented here allow describing any deviation of the final structure from symmetry. Also, the analysis of non-equilibrium structures ("failed" particles, formed under inappropriately selected process conditions) is possible using these methods (Jabłczyńska et al., 2018).

Monte Carlo methods
Monte Carlo (MC) algorithms are another group of numerical methods which found their application in numerical simulations of structured particles formation. The beginning of the development of these numerical methods had its place in the first half of 20 th century. However, they became most popular after the dissemination of computers in applications for scientific computing. They found an application not only (and not even most of all) in the simulation of particle formation but also in many branches of physics, engineering, biology and even social sciences (e.g. for modeling of the dynamics of political views in society).
In a broader sense of the term, MC methods are understood as methods that use repeated random sampling to obtain numerical results (Metropolis et al., 1953;Metropolis, 1987). The concept that stands for this method boils down to the observation that the phenomena occurring in systems with many degrees of freedom can be described with a good accuracy as random phenomena. A trivial example is the modeling of diffusion by considering the random walking of test particles. By tracking these particles, one may find out, e.g. the dependence of their mean square distance on time and thus the diffusion coefficient of the particle. The extension of this simple method is Kinetic Monte Carlo (KMC) method which is intended to simulate the time evolution of physical process (Voter et al., 2007). This evolution is realized as a sequence of "jumps" between individual states of the system. The probability of all possible "jumps" is known and thus, the "jumps" are chosen randomly. Such a sequence of states may give a picture of exemplary realization of the evolution of a process. By repeating the simulation as described above, it is possible to describe the average dynamics of the studied phenomenon and estimate possible deviations from these average dynamics resulting from the random nature of its description.
The very first applications of MC method in particle formation modeling concerned particle agglomeration. Tandon and Rosner (1999) applied a very simple MCbased scheme to find the final size distribution of agglomerates. In the mentioned work, the probability of merging particles/agglomerates has been set as proportional to the coagulation kernel (obtained from collision theory). In every MC step, the pair of particles has been randomly chosen, then a random number from 0 to 1 has been generated and compared with the probability of merging. If the random number was lower than the probability, the aggregation took place. This simple method allowed obtaining the size distribution function, which was in good agreement with the one received by means of experiments or other numerical methods. Other papers presented the application of MC algorithms for modeling the synthesis of nanoparticles from microemulsions (Tojo et al., 1997a;1997b; and the formation of the shell in core-shell nanocrystals in reverse micelle systems (Shukla and Mehra, 2006). All these early works focused on obtaining the size distribution function of the particles rather than on their internal structure.
Kinetic Monte Carlo method has also been applied for modeling the formation of structured particles. Gorshkov et al. (2014) have modeled shell formation of core-shell noble metal nanoparticles with a core made from gold and a shell-from silver or platinum. Based on their earlier work (Gorshkov et al., 2011), authors have considered the diffusion of point-like "atoms" (in general, standing for atoms, ions, or molecules) from the solution to the surface pre-formed core where they were adsorbed to form a shell. As a result of these simulations, the mean values of the parameters of particles (e.g. shell thickness) have been obtained. The numerical results appear to be in good agreement with experimental ones.
Another very popular group of MC models which have found their application for modeling of formation of NMPs are the Markov chain MC methods (Gilks et al., 1995;Andrieu et al., 2010). The clue of this conception is to obtain a chain of states of the system so these states have the desired probability distribution. For typical physical or engineering problems, this equilibrium distribution is the Boltzmann (or canonic) distribution which describes the probability of appearance of a specific state of the system as (Creutz et al., 1983): where s denotes the vector of parameters of state, E s -energy of this state, k B is the Boltzmann constant and T-temperature of the system. Generation of a Markov chain which has the canonic distribution may be done in . If a new state s 1 would not be accepted, the previous state s appears once again in the Markov chain. Then the whole procedure is repeated. Thus, it can be easily concluded that the states with a relatively high probability p B appear in Markov chain more often than those with low p B .
The works in which Markov chain MC methods were used to model the formation of particles with internal structure has appeared since the beginning of the 21 st century. The synthesis of composite latex particles of various morphology by means of two-step emulsion polymerization has been simulated using an MC algorithm (Duda and Vázquez, 2005). The authors considered a mixture of two kinds of polymers, A and B and water particles W. The interactions between these three types of individuals are described by the very simple scheme: where U ij denotes the potential energy between interacting particles, r-the distance between them, indexes i and j denote polymers A, B or water W. Finally, σ ij is a diameter of particle i if i = j or, for i ≠ j is defined as σ ij = 0.5(σ ii + σ jj )(1 + ∆ ij ). The parameter ∆ ij > 0 for every pair i and j favors the grouping of particles of the same type and thus separation of the three compounds of the mixture. To complete the description of the system, the circular domain is introduced and the energy of every particle which leaves the domain is also assigned as infinity. The MC simulations are performed in the classical canonical ensemble. As a result, these simulations lead to (dependently on the concentrations of all the compounds) to various morphologies of structured particles: core-shell, sandwich-like, raspberrylike and some other.
When investigating the final structure of structured NMPs rather than the dynamics of the formation process, MC models have appeared to be more efficient than any others what has been shown in (Chandross, 2014). In this paper, the synthesis of core-shell Cu-Ag nanoparticles has been investigated. Two different methods have been applied. The first one was MD which has been used to track every individual Cu and Ag atom during the synthesis process. The second one was the MC model. The MC method performed random particle swaps (understood as the exchange of Cu and Ag atoms at each MC step) with a small random translation below 0.2 Å. The successive states have been used to construct the Markov chain in the manner described above. The simulations show that the final state of a particle was achieved in a shorter time with the MC method than with the molecular dynamics.
While most of the applications of MC algorithms concern core-shell particle formation-this group of methods is also suitable for modeling the formation of other composite particles. In the article by Sdobnyakov et al. (2020) various morphologies of composite CuNi particles obtained by the modified solution combustion synthesis in the air have been analyzed. The MC simulations showed the patterns of neck formation for various cases of the initial arrangement of copper and nickel nanoparticles. The final structures obtained numerically were again in good agreement with XRD measurements.
MC methods have also been used for modeling the formation of porous particles from self-assembly of amphiphilic block copolymers inside an oil-in-water emulsion droplet (Zheng L. et al., 2019). The authors of this work used the MC method to obtain the final structure of such particles what allowed to construct the morphological diagram of particles as a function of the surfactant concentration and the copolymer composition (represented by the volume fraction of the hydrophilic block). They observed, dependently on simulation parameters, various morphologies of particles: non-porous, particles with closed pores and opened ones.
Among the many modifications of MC, it is worth noting the dual population balance Monte Carlo method, which allows tracking the distributions of two kinds of particles present in the system. In Skenderović et al. (2018) this method has been applied to simulate the particle formation by means of flame spray pyrolysis, where the solid particles may form from the gas phase or by means of droplet drying.
As shown, the MC models are widely used to investigate structures of NMPs synthesized in many processes. The main field of application of this group of methods is the multi-component systems. MC models are then used for the investigation of the formation of, e.g., composite or core-shell particles. When we are interested in the structure of resulting particles, the group of MC methods seems to be the most suitable tool. However, using these methods makes it difficult to track the dynamics of the synthesis and determine the structure of transitional states. This issue may be approached to a certain extent by the use of kinetic Monte Carlo methods. When we are interested in synthesis dynamics rather than in final structures, other numerical methods are worth considering.

Other approaches
The previous sections discussed the most popular numerical models for the synthesis of structured NMPs. In addition to these popular models, each of which has been described and used in at least a few scientific papers-there exist also some other models of formation of structured particles. These models were usually less versatile and were constructed ad hoc to describe a specific process that resulted in forming particles with a well-defined structure. The possibilities of generalizing such models to other cases of organized particle synthesis are relatively little. In this section, one of these less common models will be presented because of their originality, clever construction or the importance of the process they describe.
Specific methods have to be applied for the simulation of the synthesis of HPs by means of colloidal templating methods. As already mentioned, the starting point for these methods is the synthesis of CSPs followed by the removal of the core using physical phenomena such as Ostwald ripening or the Kirkendall effect. The first of these phenomena is the mass transport from smaller particles (grouped within the core) to larger ones, which make up the outer layers of the primary particle (Zeng, 2007). This process is theoretically described by the following continuity equation: Various methods of solution of Eqn. (15) are reviewed in (Baldan, 2002).
Kirkendall effect is the motion of the interface between two substances (usually metals) that occurs as a consequence of the difference in diffusion rates of the atoms of both kinds. This leads to the formation of vacancies between core and shell and shrinkage of the core and finally HPs (or yolk-shell particles). Therefore, since this phenomenon is based on diffusion processes, the Kirkendall effect models boil down to numerical solving of the diffusion equations of both components under geometrical conditions imposed by the shape and size of the originally formed CSP. From a numerical point of view, applied numerical schemes and methods are typical, e.g. finite difference methods (Gusak and Tu, 2009;Jana et al., 2013), the variational method of calculating the chemical potential of the atoms diffusing along with the core-shell boundary (Klinger et al., 2015) or the thermodynamic extremal principle which minimizes Gibbs free energy (Svoboda et al., 2009). Energy considerations have also been applied to the investigation of random shape fluctuations during the formation of HP from CSP by Kirkendall effect in room temperature (Erlebacher and Margetis, 2014). It has been (qualitatively only) shown that surface thickness fluctuations of the shell are highly probable at room temperature and may lead to exposure of the core as a result of arising of some pinholes. When core dissolves, pinholes quickly close and classic HP is formed.
For considerations of process dynamics and final structures of porous particles, the class of models which may be named "topological" has been formulated. This class is based on elementary laws of topology. Papers presenting this approach (Gradoń et al., 2004;Lee et al., 2010) investigate the formation of porous silica particles by means of spray drying of a mixture of silica and polystyrene latex (PSL) colloids above the decomposition temperature of PSL. PSL nanoparticles are much larger than silica ones and during the first stage of drying, they take such positions that the mutual distances between them are as large as possible. As long as their number in a single drop is properly selected, they are the vertices of regular polyhedron, e.g. when there are four PSL particles in a dried drop, their dots will form a regular tetrahedron. In the second stage of drying, pores (cavities) will be formed in these places due to the thermal decomposition of PSL.
Modeling of the formation of virus capsids has been developed since 1990s. First works in this subject (Zlotnick et al., 1994) applied the formalism similar to PBEs to describe the dynamics of merging of protein subunits to form icosahedral capsid. The described model assumes that the protein aggregates may join or dissociate. However, at each stage the appropriate value of the angle between the subunits is maintained-the same as the value of the doublewall angle in the icosahedra. Based on this approach, the kinetics of the concentration of all transition structures and complete capsids in time is determined. This approach to modeling the kinetics of capsid formation dominates the works in this field -sometimes KMC methods are used instead of PBEs (Kivenson and Hagan, 2010). In both cases, it is crucial to know the energy of interaction between protein subunits. For the formation of viruses with different shapes of capsid, e.g. tobacco mosaic virus, models more suited to their topology were elaborated. E.g. for the tobacco mosaic virus mentioned above, it is the so-called "zipper model" (Kraft et al., 2012). This model assumes the elongation of rod-like virus capsid by random binding of protein subunits along the helical line. In recent years works dealing with multiscale modeling of capsids formation  have appeared. Besides protein subunits merging modeling, these models also consider the protein conformation changes, what is done by means of MD-based methods.
Recently, multi-scale models and models combining several approaches described earlier have been gaining popularity. E.g. Buesser and Pratsinis (2011) proposed the simulation of gas-phase synthesis of CSPs (TiO 2 -SiO 2 ) by coupling of CFD with tracking of single particle motion. An example of the application of multiscale model may be found in (Xiao and Chen, 2014). For simulation of synthesis of two components particles in the spray drying process, the molecular-level model has been coupled with a continuum diffusion approach. The latter description gave the information about the segregation of components in the particle, while the microscopic one allowed describing the composition of single molecules at the particle surface. Similar multi-scale model, basing on CFD has been also applied for modeling of TiO 2 structured nanoparticles by means of FSP (Torabmostaedi and Zhang T., 2018).

Conclusions and perspectives
Along with the progress of research on structured particles, their synthesis and application, the development of numerical methods and models is progressing, allowing predicting the kinetics of formation and properties of these particles. This subject is still new, which is due to the fact that numerical modeling of the synthesis of structural particles, firstly, required appropriate computing power of computers, and secondly-the development of proper methods of particle characterization, enabling the verification of the results of numerical simulations. In recent years, many numerical descriptions of the synthesis of particles with a given structure have been created and are still being prepared.
The most popular are (i) methods based on the numerical solution of the differential equations of mass and energy, (ii) closely related multi-layer methods, (iii) the molecular dynamics method and other methods allowing for tracking individual atoms, groups of atoms or molecules joining a particle, (iv) methods based on population balance equations; and (v) Monte Carlo methods. The choice of the method that will be used to describe and optimize a specific process for obtaining structured particles depends mainly on the properties of the particles to be described in this process. Table 1 presents the possibilities of using the five most popular methods mentioned above to describe the synthesis of structural particles of various structures. As this table shows, molecular dynamics method and similar methods seems to be the most universal. Monte Carlo methods have an almost equally wide range of applications. It should be expected that in the near future, these methods will be most often used to simulate the NMPs synthesis processes and will be most intensively developed.
The greatest disadvantage of using these two most popular methods, i.e. the relatively long computation time and the high degree of complexity of the mathematical apparatus, has prompted and will probably continue to encourage scientists to develop other, less universal, but simpler models for the synthesis of particles with a specific structure. Models are used primarily by practitioners and experimenters to support their experiences. It is worth continuing to use these types of models where possible, being aware, however, that their use is limited to the assumptions made.

Acknowledgement
The author would like to thank Katarzyna Jabłczyńska of the Warsaw University of Technology for her valuable discussions on this topic.

Nomenclature
Abbreviations CFD computational fluid dynamics CSP core-shell particle DPD dissipative particle dynamics FSP flame spray pyrolysis HMB heat and mass balance HP hollow particle JP Janus particle KMC kinetic Monte Carlo method L-J Boltzmann constant (k B = 1.38064852 × 10 -23 J/K) L particle internal coordinate (in PBE) -size (m) L h heat of evaporation (J/kg) m mass (kg) ṁ v mass flux (kg/s) N particle number density function (m -3 ) P pressure (Pa) p B probability (-) R normalized radial coordinate (-)