Theoretical Modelling of Ion Exchange Processes in Glass: Advances and Challenges

In the last few years, some advances have been made in the theoretical modelling of ion exchange processes in glass. On the one hand, the equations that describe the evolution of the cation concentration were rewritten in a more rigorous manner. This was made into two theoretical frameworks. In the first one, the self-diffusion coefficients were assumed to be constant, whereas, in the second one, a more realistic cation behaviour was considered by taking into account the so-called mixed ion effect. Along with these equations, the boundary conditions for the usual ion exchange processes from molten salts, silver and copper films and metallic cathodes were accordingly established. On the other hand, the modelling of some ion exchange processes that have attracted a great deal of attention in recent years, including glass poling, electro-diffusion of multivalent metals and the formation/dissolution of silver nanoparticles, has been addressed. In such processes, the usual approximations that are made in ion exchange modelling are not always valid. An overview of the progress made and the remaining challenges in the modelling of these unique processes is provided at the end of this review.


Introduction
Ion exchange in glass has been used for centuries for the purposes of decoration and colouring. Glass lustre on ceramics with metallic nanoparticles from ion exchange has been known from the early Islamic culture during the 10th Century [1]. However, the scientific and industrial application of this technique dates back 60 years ago when potassium ion exchange (IE) was first applied in the chemical surface tempering of glasses [2,3]. Next, with the introduction of the concept of integrated optics in 1969 [4], ion exchange in glass was proposed as a waveguide fabrication process. Just a few years after this proposal, Izawa and Nakagome published the first work on ion exchange waveguides [5]. This kind of waveguide presents several advantages: fibre compatibility, low propagation losses and low cost. Moreover, ion exchange can be combined with other techniques, such as sol-gel, for the fabrication of passive and active (rare-earth-doped) integrated optical devices [6]. Currently, the main applications of ion exchange are in glass strengthening [7][8][9][10][11] and in the fabrication of photonic components for both guided-wave [12][13][14][15][16][17][18] and bulk optics [19][20][21].
In an IE process, cations (mostly Na + ) close to the glass surface are replaced with other monovalent cations such as K + , Li + , Rb + , Cs + , Tl + or Ag + [22]. Molten salts of such cations are common sources of dopants, although a metallic film deposited on the glass surface can be also a source of cations. The exchange takes place by a purely thermal diffusion process or it can be assisted by an electric field. The new cations can change the electrical permittivity, the stress and even the absorption of the glass, these being changes proportional to the dopant cation concentration. By selective masking of the glass surface, ion exchange can be locally prevented or allowed, giving rise to custom-made elements for specific purposes. The theoretical modelling of the IE processes ran parallel to the development of the technologies. This modelling was fundamental to the design and fabrication of most elements. For instance, in the field of integrated optics, some subjects as fibre compatibility or coupling losses depend strongly on the permittivity distribution and, hence, on the cation concentration. Therefore, the prediction of the cation concentration is of great importance to design devices based on IE technologies.
Here, we present a review of the theoretical modelling of IE processes in glass made up to now, as well as the last advances in this matter that have been reported in recent years. Ion exchange within a glass network is governed by diffusion and drift processes as a response to a concentration gradient and an electric field, respectively [23,24]. This gives rise to a concentration profile of the exchanging ions, which depends on the processing conditions of the substrate: temperature, exchange time, applied electric field, etc. The Nernst-Planck drift-diffusion equation describes this process. It establishes the proportionality among the flux density of each ion species and both the electric field and the concentration gradient. On the other hand, Poisson's equation and continuity equations for each ion must be fulfilled. This gives rise, in general, to three second-order coupled partial differential equations whose solution provides the evolution of both the cation concentrations and the electric potential. In the derivation of these equations, the charge neutrality approximation is usually assumed, which allows for some important simplifications without a relevant loss of accuracy in the most common cases [25,26]. All of these subjects are addressed in Section 2, where we present the basic model of ion exchange. This model assumes that the self-diffusion coefficients of the exchanged cations are constant for a given temperature. However, experimental measurements [27] showed that these coefficients depend on the cation concentrations. Therefore, this basic model was generalized to non-ideal cation behaviour and concentration-dependent self-diffusion coefficients by considering the so-called mixed ion effect [28,29]. Later, both models were rigorously generalised, and in the derivation of the equations, Faraday's law was considered instead of Ohm's law [30]. This leads to a non-standard Laplace equation for an effective electric potential. Accordingly, the boundary conditions for the most common IE processes were established. These conditions, together with the aforementioned partial differential equations, complete the theoretical modelling of the IE problem. On the other hand, some IE processes have attracted a great deal of attention in the last few years due to their remarkable applications. Among these, we must highlight: glass poling, electro-diffusion of multivalent metals and the formation/dissolution of silver nanoparticles. Their modelling has only been partially done so far [31][32][33], because the usual theoretical assumptions (mainly charge neutrality approximation and ideal cation behaviour) are not always valid in such processes. In Section 4, we give an overview of the progress made and the challenges still to be faced on this matter.

Basic Model
The simplest problem of ion exchange arise when two species of monovalent cations (A and B) exchange with each other in the same glass region at a given temperature T. Let us consider an infinite one-dimensional medium (x being the spatial coordinate) with a homogeneous concentration C 0 of fixed anions (typically -Si-O − radicals) and nonhomogeneous and variable cation concentrations C A (x, t) and C B (x, t), t being the time. Moreover, we considered that each cation is initially near its anion, that is they are paired, so C A (x, 0) + C B (x, 0) = C 0 . Consequently, space charge density is initially cancelled: that is local charge neutrality is met. Our goal was to obtain the evolution of cation concentrations given these initial conditions.
where J i is the flux density and D i the diffusion coefficient of each cation. However, the two interdiffusing cations have usually different diffusion coefficients and, therefore, different mobilities, which produce charge imbalances. These imbalances generate a strong internal electric field (E), which tends to balance the charges and restore charge neutrality (Equation (1)). Therefore, Fick's equation is no longer valid, and a drift term must be added on its right-hand side to take into account the effect of this electric field on the cation motion. This leads to the Nernst-Planck equation [34]: where e is the proton charge, T the absolute temperature and k Boltzmann's constant. Note that some authors included, in the drift term of this equation, the Haven ratio. However, as we will see below, this parameter should not be incorporated into the model as a general rule. The occurrence of the above-mentioned electric field can be seen from a quantitative point of view by calculating the total flux density: which depends on the electric field through mobility u: However, J 0 cancels in the current problem, that is, because there is no external field applied that generates a net current. This means that the imbalance of the diffusion of the two cations (D A ∇C A + D B ∇C B ) is compensated by the internal field through the drift term uC 0 E. Note that Equation (4) is a generalization of Ohm's law. On the other hand, as long as there is no creation or destruction of cations from/to a metallic state, the continuity equation must be fulfilled: Now, by combining this equation and Equation (3) for, for instance, i=A and substituting, in the resulting equation, the electric field from Equations (4) and (6), we obtained the differential equation that gives the concentration evolution [23] of cation A: where the charge neutrality (Equation (1)) was applied and a normalized concentration (c A = C A /C 0 ) was used, and we defined the following interdiffusion coefficient: where α = 1 − D A /D B . The dependence of u andD on c K predicted by this basic model, for K + /Na + IE, can be seen in Figure 1a.
Self-diffusion and interdiffusion coefficients, as well as mobility, calculated from experimental data obtained by radiative tracers [27], as a function of the normalized concentration for K + /Na + IE. The Haven ratio was ignored. (a) Basic model-Equations (5) and (9)-which assumes that the self-diffusion coefficients remain constant with the cation mole fraction. (b) The same functions taking into account the MIE. Quadratic polynomials in c K were fitted to the logarithms of the experimental self-diffusion coefficients and then used to calculate the rest of the functions through the definitions (35) and (34). Note the difference betweenD and D mob , which are the expected interdiffusion coefficients when the interaction among cations or the ideal mixture is assumed, respectively.
A more complex problem is the IE assisted by an external electric field. In such a case, Equation (6) is not met. Therefore, additional equations are necessary to calculate the total flux density (J 0 ) or, alternatively, E, which now includes the external field. As for E, it must fulfil Poisson's equation and Faraday's law of induction: where is the glass electrical permittivity and B is the magnetic field, which will be assumed as time independent since the total current changes very slowly. On the other hand, we also assumed the charge neutrality approximation (Equation (1)). This cannot be done in general, due to the existence of the aforementioned external field; however, in most cases, this is a very good approximation (see the next subsection). Doing this, the addition of continuity Equation (7) leads to: Now, by finding the electric field in Equation (4): taking into account Equation (12) and doing the same steps as before, the following electrodiffusion equation is obtained: which is an extension of Equation (8). On the other hand, Equation (13) can be inserted into (11), leading to: where we used Equation (9). This equation, Equation (12), and the boundary conditions, which will be presented later, determinate the flux density J 0 . From this flux density, Equation (14) will provide the evolution of the concentration of cations. Finally, once J 0 and c A are known, the charge neutrality approximation can be checked. This will be analysed in the following subsection. Alternatively, the last four equations can be expressed in terms of a scalar function. Indeed, an irrotational vector field is the gradient of a potential function φ, so: where the factor −e/(kT) was included in order for φ to have the same units as the electric potential. From this equation and Equation (12), we obtained a non-standard Laplace equation: which is more convenient for resolution purposes than vector Equations (12) and (15). Note that this effective potential φ is not the electric potential V, whose minus gradient is the electric field E given by Equation (13). However, a relationship between them can be obtained [30] by inserting Equation (16) into Equation (13), that is: Under typical IE conditions, the difference between V and φ is no greater than a few tenths of a volt, which is negligible compared to the usual voltages (20-100 V) used in field-assisted IE. Despite this, this difference must not be ignored in the modelling, as significant errors in the calculation of the electric field could be made. Indeed, although both potentials are similar, their gradient is not always.
On the other hand, it is worth mentioning that Equations (12) and (15) are trivially solved in the one-dimensional case, that is J 0 is constant. This constant will be established from the experimental setup. If a constant current source is used, this value is set directly.
Otherwise, when the external field is generated by a constant voltage source, the value of J 0 can be calculated from the voltage applied to the sample and Equations (14) and (16)-(18).

Charge Neutrality Approximation in Field-Assisted Ion Exchange
In the previous subsection, the charge neutrality was assumed in the field-assisted IE problem, as well as in the thermal-only case. However, this charge neutrality is not always fulfilled, especially when strong external electric fields are used and/or very high concentration gradients exist. In fact, some authors have modelled the silver concentration in channel waveguides by considering explicitly the space charge distribution [35], albeit at the cost of including other assumptions.
An estimation of the validity of the charge neutrality approximation can be made from Equation (13), which can be expressed as: where we used the charge neutrality approximation and the definition of α. Now, we took divergences in this equation in order to compare it with Equation (10) and estimate the error made by this approximation. This is an iterative procedure. First, we used this approximated expression for the electric field. Next, we substituted it into Poisson's Equation (10) to obtain a more accurate value for C A + C B , which was introduced in the previous equation, and so on. However, a unique iteration will be enough to obtain an order of magnitude of the charge density [36]. Therefore, if we use Equation (12) and assume that does not depend on C A and C B , we obtained: Now, from this equation, we can obtain some conditions for the validity of the charge neutrality approximation by comparing each term on the right-hand side of this equation with the second term on the left-hand side. Therefore, for the second term on the right-hand side, this comparison provides: for a BK7 glass with a density of 2.4 g/cm 3 and 8.4% by weight of Na 2 O, which gives C 0 3.9 × 10 27 m −3 ; likewise, typical values for the temperature, T = 400 • C, and for the relative permittivity, r = 10, were considered. As for the ion exchange, we chose for this assessment a Na + /K + IE in a soda-lime glass with a diffusion coefficient ratio D K /D Na 2.5 × 10 −2 [37], which leads to α = 0.975. Therefore, ∇c A must fulfil the following condition: In the worst case, c A = 1. This condition means a change, in the normalized concentration, from one to zero-point-nine at a distance much greater than 0.36 nm, which in practice is met in almost any IE process. Indeed, most IE processes use glasses with a higher Na 2 O content or their α values are less than the one of Na + /K + IE. Therefore, this condition is even less restrictive in those cases.
On the other hand, the third term on the right-hand side of Equation (20) is small under the same conditions as the second one. Finally, the first one is small if the second one is and, furthermore, if: where we took into account [37] that D B 6 × 10 −13 m 2 /s, for T = 400 • C. That is, in the worst case (c A = 1), the current density must fulfil: This condition is ensured in practice because the highest current densities that have been used in field-assisted IE processes until now are much lower than this value. In fact, Joule heating should be taken into account for currents higher than a few mA/cm 2 , due to the strong dependence of diffusion coefficients on the temperature [38,39]. Therefore, the charge neutrality approximation is valid in most IE process, even for Na + /Rb + and Na + /Cs + ion exchanges, which give rise to very steep concentration profiles. In these kinds of profiles, although ∇c A may be very high, the denominator of Equations (21) and (23) is not close to zero, even for molar fractions close to one. Indeed, we assumed that α = 1 − D A /D B is constant, but actually, the diffusion coefficients depend on the molar fractions, so that if a cation has a low concentration, it also has a lower diffusion coefficient. In other words, if c A 1, then α < 0, and 1 − αc A is never close to zero. The basic model that we present in this section cannot explain these dependencies on the concentration of the diffusion coefficients, and therefore, this cannot be used to describe the whole range of molar fractions.

Boundary Conditions
When an external medium is in contact with the glass surface, each type of cation present in the glass and/or in the medium can either cross or not cross the glass/medium interface. This depends on the nature of such a medium. If a cation crosses the interface, a thermodynamic equilibrium is assumed between the crossing cations at both sides. Otherwise, the normal component of its flux density cancels. Therefore, the different boundary conditions for the Equations (14), (16) and (17) are obtained for a mask/air, a silver or copper film and a molten salt mixture, when none, one or both cations, respectively, cross the interface. Moreover, when the external medium is a conductor (metallic film or fused salt), its electric potential can also be chosen and directly affects the boundary condition of the potential φ through Equation (18).
In Tables 1 and 2, we summarize the boundary conditions for IE processes from molten salt mixtures and films, as well as under a mask or in air. The equations that govern these processes are also shown. Note that boundary conditions for other common IE processes as annealing or secondary ion exchanges are included in these cases. Indeed, boundary conditions for annealing processes are the same as for "mask/air", and the ones for secondary ion exchanges are included in the "salt" case, just using a different constant "C", which depends on the dopant concentration in the salt mixture and temperature through an equilibrium equation [40][41][42][43]. Moreover, the final concentration of the first process is the initial condition for the second one.
We show, in these tables, the two alternative forms presented in Section 2.1. In Table 1, the problem is formalized in terms of the normalized concentration of dopant cations c A and the total flux density J 0 . This is a more straightforward form, which makes it easier to understand the physical problem of ion exchange. However, for resolution purposes, the form given in Table 2 is simpler because it is written in terms of c A and a scalar function φ, which can be seen as an effective electric potential [30].
The basic model of this section implicitly assumes an ideal cation behaviour in the glass. Therefore, for the sake of consistency, the same assumption was made in the electrochemical potentials to obtain the theoretical equations for the boundary conditions of Tables 1 and 2. However, the corresponding condition for c A in a glass-salt interface, as a function of the salt composition, must contain non-ideal terms on the glass side, as we will show in the next section. Owing to the molten salt being homogeneous, c A | S is constant along the surface. Hence, a mixed theory can be a pragmatic approach; similarly, experimental values for c A | S can be used to feed the numerical algorithms that model the ion exchange. On the other hand, the boundary conditions for φ are sometimes irrelevant (thermal diffusion) or they can be approximated to those of V (diffusion assisted with strong fields). As a result, they were often ignored in models and experiments in the context of waveguide fabrication. Table 1. Complete formulation of the IE problem (equations and boundary conditions), for the most common IE processes, in terms of the normalized concentration of dopant cations c A and the total flux density J 0 ;ê S is a unit vector perpendicular to the sample surface, and C is a constant, which mainly depends on the dopant concentration in the salt.

Some Particular Solutions
Two noticeable one-dimensional solutions are obtained for both specific thermal and field-assisted IE processes. In the simplest thermal IE, a glass sheet is immersed in a molten salt, which contains foreign cations. The boundary conditions for c A and φ are constant over time and along the glass surface, then J 0 = 0. The concentration profile of the dopants scales in depth proportionally to the square root of the diffusion time, but it retains its shape, that is c( , being the origin of the coordinates at the glass surface [44]. The shape of f depends on c A | S and α. If either of them is small, the diffusion Equation (14) becomes linear and f tends to a complementary error function (erfc), which has an inflection point at the glass surface and decreases monotonically to zero inside the glass. Otherwise, the diffusion is faster near the surface, where the dopant concentration is high. This causes f to present a bump near the surface, being the inflection point displaced into the glass. f is often approximated by a Gaussian function, although other functions have been proposed [45][46][47].
When a voltage is applied between both surfaces of the glass sheet, a current normal to these surfaces appears. If J 0 is kept constant, a stationary solution of Equation (14) exists when the slow cations invade a region containing a higher concentration of fast cations [24]. In that case, the mole fraction profile of the slow cations is a step-like function that moves at a constant speed while maintaining its shape. Therefore, two regions with different concentrations are formed with an abrupt transition between them. The higher the step is, the sharper the front and the greater its velocity because the mobility of slow cations increases with its concentration. The stability results from a competition between diffusion and the non-linearity of the drift term. The former widens the profile, whereas the latter sharpens it because the rear region of the profile moves faster. On the other hand, if V is kept constant, J 0 decreases slowly over time, since the resistivity of the sample increases as the doped region becomes thicker. This effect may already be appreciable for doped regions a few microns deep, because the mobility ratio can reach several orders of magnitude [38,48]. In [38], the applied voltage V had to be corrected by 0.93 V to accurately describe the experimental depths because these models did not include the potential φ. Note that a combination of the boundary conditions of φ at both glass sides, as well as the difference between φ and V given by Equation (18) explain such a potential. Finally, it is worth noting that a stable profile is only formed if the slow cation chases the fast one. Otherwise, the profile extends indefinitely.

The Mixed Ion Effect
Constant diffusion coefficients were assumed in Section 2. However, they depend on both the temperature T and the molar fraction c A . In monoalkaline glasses, the temperature dependence follows the Arrhenius law: This is due to each cation needing to surmount a potential barrier of height equal to Q to move to another potential well. The diffusion coefficient is proportional to the number of cations that overcome this energy in a given instant. This number follows a Maxwell-Boltzmann statistic, which explains the above exponential law. According to the Stuart-Anderson model [49], the potential barrier has both an electrostatic component, corresponding to the energy necessary to separate the cation from its anion, and a mechanical component that describes the glass network distortion necessary for the cation to break through another potential well. Therefore, when two species of cations are present in the glass, it is natural that each one has its own diffusion coefficient. However, surprisingly, the D 0 and Q values of each species (D 0A , Q A and D 0B , Q B ) largely depend on the cation mole fraction, in such a way that each diffusion coefficient is reduced by up to several orders of magnitude as its respective mole fraction is approaching zero [27,50,51]. This reduction is mainly due to an increase of the activation energy of minority cations; furthermore, it is stronger than the difference between the diffusion coefficients of both species in their respective monoalkaline glasses (D AA ≡ D A | c A =1 and D BB ≡ D B | c A =0 ). Consequently, D A and D B are equal for some intermediate mole fractions. We will show that this leads to a maximum value of the interdiffusion coefficient and a minimum value of the direct current (DC) conductivity for this intermediate mole fraction or its neighbourhood. Therefore, the DC conductivity shows an excess of activation energy with respect to the monoalkali glasses (this excess disappears with frequency in AC conductivity). This phenomenon, among others, is included in the so-called double-alkali effect, mixed alkali effect or, later, mixed mobile ion effect or mixed ion effect (MIE), since silver, thallium or copper cations behave similarly to alkali ones in glass [37,[52][53][54].

Brief Review of Theories
The origin of the MIE has been debated for a long time, but no theory has been universally accepted yet. Some theories focus on an interaction among neighbouring cations in such a way that mixed pairs are assumed to be more energetically stable than pairs of cations of the same species [15,[55][56][57][58][59]. This assumption is supported by several experimental achievements. Namely, the interdiffusion coefficient presents a thermodynamic term for alkali IE [60,61] and for Ag + /Na + exchange [62]; we will show that this term is missing if cations behave as an ideal gas. Furthermore, the surface mole fraction of cations in ion exchanged glasses from molten salts (the salt boundary condition) must be explained on the assumption that a non-ideal cation behaviour both for double-alkali exchange [40,41] and for Ag + /Na + exchange [42,43]. On the other hand, mixing enthalpy experiments do not have a clear interpretation. A negative mixing enthalpy (net attraction among dissimilar cations) was found [63][64][65], which correlates linearly with the excess of activation energy for DC conduction. However, the former is about 20 times weaker than the latter. Besides, mixed pairs of cations would be expected to be much more likely than pairs of the same species in the presence of interaction, but nuclear magnetic resonance experiments [66][67][68], as well as neutron and X-ray diffraction [69] show that they are rather randomly distributed. Consequently, other authors attributed the MIE to relaxation processes in the glass structure [70][71][72][73]. In particular, each cation is assumed to modify its neighbourhood after a relaxation time to achieve an energetically favourable site. Therefore, in a monoalkaline glass, all the sites are of the same type. Once a foreign cation enters this glass and modifies a site, it is difficult for it to diffuse because all accessible sites are the wrong type. Even if the cation gains access to one of them, most likely, it will return before the new site relaxes. Moreover, diffusion is assumed to occur through conduction pathways, so a foreign cation can block several indigenous ones. This explanation is compatible with the ideal behaviour of cations, that is with a random distribution.
The theory that we assumed is relevant because, depending on it, the resulting diffusion equation is slightly different, as we will show below.

The Cation Flux Density
Let us consider the electrochemical potential (μ i ) of each cation species that has a chemical term depending on the thermodynamic activity (a i ) of that species and an electric term proportional to the cation charge (e) and the electric potential (V): The flux density of each cation species is proportional to both the gradient of its electrochemical potential and its concentration. This leads to the Nernst-Planck equation: where the g i 's are the thermodynamic factors: As mentioned above, some authors include the Haven ratio in the drift term of Equation (27). However, the Haven ratio should only be used to obtain D i from the experimental values of the diffusion coefficient of radioactive tracers: D * i [74]. The difference between them arises, for example, if the diffusion mechanism is the indirect interstitial one. In this case, a cation located in one interstice replaces a nearby regular site cation, which jumps to another interstitial site. Consequently, the total mass (or charge) displacement described by D i is different from that of a single cation, which is measured by the tracer. Moreover, the tracer is also affected by the thermodynamic factor, then a comparison between the mobility of a species and the corresponding tracer diffusion coefficient can result in an apparent Haven ratio [61].
If the MIE is fully due to the relaxations of the glass structure, the cation behaviour being ideal, the activity will be equal to the mole fraction, so g i = 1. On the contrary, if the cation interaction is the only thing responsible for the MIE, we will obtain D i = D ii γ i , γ i being the thermodynamic activity coefficient (a i = γ i c i ). Let us deduce that result. We divided all cations A into two sets, the ones that are hopping from one site to another in a given instant (A ↑ ) and the rest of them, which are fixed (A ↓ ). Although the former are much scarcer, both sets are in thermal equilibrium in any given small region in the glass; therefore:μ Now, we made two assumptions. First, we supposed that mobile cations behave ideally with respect to each other due their low concentration (C A ↑ C A ↓ ): Note that this could fail in the case of cooperative movement, that is, when two or more nearby cations change their site simultaneously. Second, we assumed that the reference potential of mobile cations (µ 0 A ↑ ) is independent of c A ; therefore, ∇µ 0 A ↑ = 0. In addition to using Equation (27), the cation flux density can also be calculated from mobile cations, being proportional to both the gradient of their electrochemical potential and their concentration: where D A ↑ is a temperature-dependent multiplicative coefficient and the factor 1/kT was introduced in order for D A ↑ to have units of a diffusion coefficient. By combining Equations (26), (28) and (29), we can find C A ↑ as: and replace it in Equation (30). The resulting flux density is: which agrees with Equation (27) by identifying: Obviously, an identical derivation can be done for B cations to obtain D B = D BB γ B . Surprisingly, it was not necessary to make any assumptions about the particular dependence of γ i on the mole fraction. In short, under these cation interaction assumptions, the dependence of the diffusion coefficient of each species on the mole fraction was directly related to its thermodynamic activity coefficient. For the equations to remain valid, regardless of the MIE explanation, we will continue the derivation from Equation (27), without any assumptions on the thermodynamic term g i , which can be done in the last step.

Generalized Equations and Boundary Conditions
By following the same procedure as in the previous section, we replaced E with J 0 in the expression of J A , and then, we applied the continuity condition to obtain: where the whole flux density J 0 can be obtained from the scalar potential φ as: and φ in turn satisfies the following non-standard Laplace equation: In view of Equation (31), we can redefine the interdiffusion coefficient as: It can be split into a mobility term D mob and a thermodynamic term (g B c A + g A c B ), which is not in the definition (9) of the basic model. If cations behave ideally, the latter becomes equal to one, and the interdiffusion coefficient is reduced to the mobility term. If cation interactions are relevant, the thermodynamic term enhances the maximum ofD(c A ) for intermediate values of c A , as can be seen in Figure 1b. Similarly, Equation (32) shows that the mobility is: which, in fact, is the same as Equation (5), but now, the D i 's are mole fraction dependent. This results in a minimum conductivity value for a mole fraction close to that at which the interdiffusion coefficient reaches its maximum. In contrast, the basic model leads to monotonic u andD functions. In order to impose the boundary conditions on φ, we need to relate it with the electric potential v. This can be done by the procedure followed in the above section, but the resulting expression is not so simple: If the MIE is caused by cation interactions (D i c i = D ii a i ), this expression can be integrated for any particular dependency of the activities on the mole fraction: On the contrary, if cations behave ideally and relaxation processes are responsible for the MIE, then Equation (36) becomes: and the integration can only be performed once the dependencies of the diffusion coefficients on the mole fraction are known. The starting point to impose the boundary conditions are the same as in Section 2, that is equal electrochemical potentials at both sides of the glass surface for cation species that can cross it and zero flux through it otherwise. However, the particular functions for the electrochemical potentials and flux densities must be chosen from the model of the glass behaviour.

Changes in the Solutions with Respect to the Basic Model
The main difference, with respect to the basic model, which is observed after the numerical solution of Equation (31), comes from the presence of a maximum in the interdiffusion coefficient for intermediate molar fractions (Figure 1b). Therefore, new qualitative shapes of the profiles are seen when almost all indigenous cations are replaced by the foreign ones.
In Figure 2, we show a simulation of molar fraction profiles of potassium cations in a K + /Na + thermal exchange for several boundary conditions ( c K | S ). The profile for the lowest value of c K | S was similar to an erfc function. For intermediate values (0.4 and 0.6), the profiles showed a bump between the glass surface and the tail, near which there was an inflection point. Both characteristics were similar to those predicted by the basic model. Nevertheless, a second inflection point (or equivalently, a high slope at the surface) appeared for the highest values of the boundary condition. Because the interdiffusion was, in the present model, lower near the surface, the slope of the profile must increase to provide cations at a sufficient rate for the intermediate regions, where they diffuse faster. An approximate analytical profile was proposed to describe all of these cases [75]. In this work, the authors checked the quality of that profile through the measurement of the effective indices of ion exchange waveguides. They obtained an average deviation of 3 × 10 −4 between measured and calculated effective indices, which corresponded to a deviation of ≈0.3% between concentration profiles. Similar results have been obtained by numerically solving the diffusion equation that governs the IE process [22]. Moreover, this author even monitored such a process as a function of the diffusion time.  A simulation for field-assisted IE can be seen in Figure 3. For the lowest boundary conditions ( c K | S = 0.1-0.3), the exchange was dominated by the diffusive term, and the stable profile was not formed yet. For intermediate values of c K | S (0.4-0.8), the stable profile was clearly formed, since the rear region of the profile was flat. Moreover, we can see that the velocity of the front and its slope increased with the boundary condition. Again, these characteristics agreed with the prediction of the basic model. Finally, for c K | S = 0.9-1, a new effect appeared. The speed and the slope of the front saturated, while its rear part was no longer flat because the potassium cations were faster here than the sodium ones. Both regimes can be observed in experimental studies [76,77].

Future Challenges
Some of the previous assumptions in the modelling of ion exchange in glass are not always valid. For example, the interdiffusion coefficient is reduced during the ion exchange of dopants, which generates compressive stress. This is probably due to a constriction of the interstitial sites [78][79][80]. Another issue is the generalization of the charge neutrality approximation to processes involving the simultaneous diffusion of three species (for example sodium, silver and potassium). Some of such processes were proposed in the past to improve the shape of narrow-channel waveguides [81,82] or, more recently, to obtain both antimicrobial and strengthening properties [10]. In [82], a simplified model was presented by considering that potassium cations are immobile. Unfortunately, little progress has been made in all of these issues in the last few years. Note that modelling of these processes in the MIE framework would require considerable experimental studies to obtain the interdiffusion coefficients. Nevertheless, there are some processes that also require more complete models, but that have recently attracted a great deal of attention due to their remarkable applications. Among them, we highlight: glass poling, electro-diffusion of multivalent metals and the formation/dissolution of silver nanoparticles.

Glass Poling
Glass poling is the distribution of the electric charges of the glass. This may result in a permanent electric field inside the glass, which can give rise to relevant non-linear effects on light propagation. Applications of glass poling include grating fabrication [83], second-harmonic generation [84] or the fabrication of optical waveguides with profiled electrodes [85]. Glass poling was initially applied to silica glass [86], which contains residual amounts of cations (C 0 ∼ 10 23 m −3 ). It is realized by subjecting the sample to a strong potential difference U (typically, a few kV) at a temperature of about 250-300 • C and, usually, not allowing charged species to enter the glass ("blocking anode") [87]. Therefore, the field rearranges the cations until a new equilibrium is achieved when the original field inside the glass is cancelled by charges located near the glass surface. That is, the applied field creates a layer, under the anode, a few micrometers thick and depleted of cations. Obviously, charge neutrality is not fulfilled here. Instead, the one-dimensional form of Equations (3) and (10) of the basic model shows that the layer thickness d, in equilibrium, is given by: This layer is charged and has a high electrical resistivity since its anions are fixed to the glass network. Furthermore, because the field is cancelled in the glass bulk, the applied voltage drops across the depleted layer. Therefore, a very strong field arises, which is independent of the thickness of the substrate. Besides, potentials above ∼1 kV generate structural changes in the layer [88]. Then, if the sample is sharply cooled to room temperature, the distribution of charges becomes frozen. Electric fields inside the silica glass can reach ∼10 7 V/m, which induces a non-linear coefficient χ(2) ∼ 1 pm/V [89]. However, this process is idealized. In practice, other cations are often involved due to the high applied voltages. One of the possible cations is hydronium (H 3 O + ) from air water vapour that forms a naturally hydrated layer in the glass. In this situation, called poling with a non-blocking anode, Equation (37) is no longer valid. Instead, a stable profile, such as described in Section 2.4, is formed because hydronium moves slower than sodium. In this case, a partially depleted region is still formed before the electric field increases enough to move the hydronium cations significantly [90]. The charge and the electric field are expected to be at their greatest values when the hydronium layer begins to form. Subsequently, a part of the applied voltage U drops across the hydronium layer [91]. In glasses with high alkali content, the situation is somewhat different, since air cannot provide hydronium at a sufficient rate. Moreover, as C 0 is large, d is small (see Equation (37)), and the electric field can become high enough to move other cations such as Ca +2 or Mg +2 [31,92]. Another situation of interest occurs when a BK7 glass, which contains sodium and potassium, is poled. Initially, only sodium cations are removed from the glass surface, as they move faster. Then, the field increases in the charged layer until potassium cations start to move and accumulates behind the sodium ones, filling the empty sites that the latter have just left. Simultaneously, a fully depleted surface layer is formed [92]. Poling of double-alkali glasses, with a non-blocking anode, was also simulated with similar conclusions [93]. One interesting result, which is obtained in high-alkali glasses, is the formation of a waveguide next to the depleted layer, this having non-linear properties. In all the previously cited works about the simulation of glass poling, the authors assumed diffusion coefficients constant, which is only an approximation. The diffusion coefficient in monoalkali glasses depends notably on the alkali content. Therefore, it is expected that the diffusion coefficient in the depleted layer is also different from bulk glass. Therefore, more experimental research on glass poling is still necessary to make clear, on the one hand, the conditions under which the charge neutrality approximation is valid and, on the other hand, whether the MIE or other effects are relevant enough to be considered in the modelling.

Electro-Diffusion of Multivalent Metals
Monovalent cations as Ag + /Na + and K + /Na + are by far the most used cations in IE processes in glass. However, some successful attempts at doping glass with multivalent cations have also been realized, mainly with transition metals, but also with rare-earths. For example, Gonella et al. obtained Co 2+ , Au 3+ and Cr 3+ diffusion profiles in silicate glasses [94,95]; and Cattaruzza et al. introduced Er 3+ cations inside soda-lime glasses [96]. These authors used field-assisted configurations with electric fields up to 400 V/mm, as thermal-assisted IE is not very effective with multivalent ion species. This allowed for penetrations of cations of up to ∼1 µm from deposited films. The main utility of glass doping with multivalent transition metals arises from the possibility of converting the glass into an active media. Therefore, optical amplification or waveguide lasers have been demonstrated in erbium-doped glasses [97], and Cr 4+ -doped materials have been proposed as both laser gain materials and saturable absorbers for passive Q-switching in IR lasers [98].
The ion exchange with multivalent cations gives rise to concentration-dependent structural changes in the glass matrix. This is due to the local coordination rearrangements, which are produced at the ion sites [99]. In addition, the amount of metal that penetrates into the glass matrix, as well as the shape of the concentration profiles depend strongly on the process parameters. Therefore, the model presented in Section 3, which includes the MIE, should be used in the modelling of such processes. Likewise, when the multivalent cations are introduced into the glass, a region depleted of cations, similar to that observed in poled glasses, is formed. Therefore, the approaches used in the modelling of glass poling could be applied to the present processes. However, little work has been done in these directions, and a comprehensive model of electro-diffusion of multivalent metals is still lacking.

Formation/Dissolution of Silver Nanoparticles
Metal nanocluster formation following ion exchange can arise if both suitable dopant cations and post-exchange processing techniques are used. Among all dopants, silver stands out for its great diffusivity in glass and its high tendency to form metal clusters [33]. As for the post-exchange techniques, heat treatments (annealing) under different atmospheres [100,101] and (pulsed) laser irradiation [102] are commonly applied. The formation of metal structures can be regarded as a drawback and an advantage. For example, in light waveguiding, metal inhomogeneities must be avoided because they cause a high light absorption. On the contrary, some investigations have found these structures useful for sensing applications, through the surface enhancement of Raman scattering (SERS) [103], for fabricating photonic crystals [104] or, in general, for plasmonic optics [105,106].
Modelling of the formation and dissolution of silver nanoparticles is a very challenging task because the two complex processes compete dynamically. On the one hand, annealing or irradiation promotes silver aggregation and the formation of nanometre-sized clusters while, at the same time, the diffusion and dissolution of silver cations are produced. Therefore, at least three species must be included in the modelling: both exchanging cations and metallic silver. Moreover, this process is often done in gaseous atmospheres, which provide other cations that diffuse into the glass, increasing the number of species to be included in the model. On the other hand, this kind of process is strongly dependent on the cation concentration, so that the MIE should be considered. Finally, not enough experimental studies have been realized until now to characterize the great variety of mechanisms involved in these processes. Some theoretical descriptions of nanoparticle formation under specific post-exchange processing techniques (purely thermal annealing [107] and annealing in hydrogen atmosphere [108]) have been given. Both authors assumed the charge neutrality approximation that, under these processing techniques, is probably fulfilled. However, they did not assess the introduction of the MIE in their models and if that would lead to a better description of the cluster formation. Considerable work, both theoretical and experimental, must be done before having a complete theoretical model of this promising technique.

Conclusions
The improvements made in the last few years in the theoretical modelling of ion exchange in glass have contributed to a better understanding of the physical and chemical mechanisms involved in this process. The basic model that has been used for decades only allows for the accurate description of ion exchange problems where the dopant concentration is low. If higher concentrations are present, a model based on the mixed alkali effect must be considered. However, some aspects of this model are still open, because the physics and chemistry involved are not fully understood. On the other hand, a significant progress has been made in the modelling of some IE processes that have attracted much attention in recent years. These processes include glass poling, electrodiffusion of multivalent metals, and the formation/dissolution of silver nanoparticles. All of them have remarkable applications, both scientific, as well as technological.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

DC
Direct current IE Ion exchange MIE Mixed ion effect