Terahertz time-domain spectroscopy for powder compact porosity and pore shape measurements: An error analysis of the anisotropic bruggeman model

Terahertz time-domain spectroscopy (THz-TDS) is a novel technique which has been applied for pore structure analysis and porosity measurements. For this, mainly the anisotropic Bruggeman (AB-EMA) model is applied to correlate the effective refractive index (neff) of a tablet and the porosity as well as to evaluate the pore shape based on the depolarisation factor L. This paper investigates possible error sources of the AB-EMA for THz-TDS based tablet analysis. The effect of absorption and tablet anisotropy – changes of pore shape with porosity and density distribution – have been investigated. The results suggest that high tablet absorption has a negligible effect on the accuracy of the AB-EMA. In regards of tablet anisotropy the accuracy of the porosity determination is not impaired significantly. However, density distribution and variations in the pore shape with porosity resulted in an unreliable extraction of the tablet pore shape. As an extension of the AB-EMA a new concept was introduced to convert the model into bounds for L. This new approach was found useful to investigate tablet pore shape but also the applicability of the AB-EMA for an unknown set of data.


Introduction
Tablets are considered the most convenient way to deliver active pharmaceutical ingredients (APIs) to a patient given their ability to achieve accurate dosing, long shelf life and cost-effective production (Desai et al., 2016). Porosity and pore anisotropy are recognised as important parameters that influence tablet disintegration and API dissolution (Markl et al., 2018a). Important precursors of disintegration such as swelling and liquid imbibition have been found to be dependent on total tablet porosity (Yassin et al., 2015a(Yassin et al., , 2015b. However, it is also well known that total porosity alone is not sufficient to fully describe the liquid imbibition process (Berg, 2014) and that pore anisotropy including tortuosity and connectivity need to be considered as well amongst other factors (Markl et al., 2018b).
Terahertz time-domain spectroscopy (THz-TDS) has been suggested as a fast and non-destructive technique for porosity and pore shape measurements (Bawuah et al., 2016b;Bawuah and Peiponen, 2016;Bawuah et al., 2016c). The use of non-ionising radiation and the advantage of fast (sub-second) measurement acquisition makes it ideal for in-line or at-line control applications (Bawuah et al., 2020). Furthermore, most pharmaceutical excipients are at least semitransparent at terahertz frequencies allowing for measurements in transmission mode (Zeitler et al., 2007). Due to the nature of the measurements being performed in the time-domain, the refractive index can be extracted from the measurements without relying on more complicated models such as the Kramers-Kronig relations, which may require further interpretation or constraints applied to succeed (Bawuah et al., 2014).
The refractive index at terahertz frequencies has been found to correlate with the porosity (Bawuah et al., 2014(Bawuah et al., , 2016bJuuti et al., 2009). In addition to being able to measure total porosity, by using different physical models, it is possible to gain further insight into the pore structure of tablets. Such advanced analyses of the terahertz data is implemented using models from effective medium theory. The anisotropic Bruggeman model (AB-EMA) is the most commonly used approach and this method has been demonstrated to allow for the extraction of pore shape information in addition to the tablet porosity (Bawuah et al., 2020).
The focus of this work will be to provide an error analysis of the AB-EMA model for THz-TDS based porosity and pore shape analysis. It will critically examine the impact and sensitivity of the following assumptions that are required when applying the model to experimental data on the accuracy of the results: I) The AB-EMA utilises the complex dielectric permittivity (ε) which is defined by ε =ñ 2 = (n − iκ) 2 , where ñ is the complex refractive index which is composed of the refractive index, n, and κ, the extinction coefficient. Thus far, all studies published in the pharmaceutical literature for using AB-EMA for pore analysis only consider n and the effect due to absorption from the tablet matrix is approximated to be negligible. II) The AB-EMA method considers pores of different porosities to have the same shape, which can be described by the so-called depolarisation factor, L. However, the tablet compaction process is inherently anisotropic given the tool and process geometry. It is also well-known that increasing compaction pressure (i.e. decreasing porosity) can result in higher tablet anisotropy (Porion et al., 2010). The elastic and plastic deformation processes themselves are strongly dependent on the magnitude and direction of applied pressure (Lubarda, 2016). Since pores in tablets are formed by the void space between particles, deformation of one is considered to affect the form of the other. III) In the most common implementation, the THz-TDS transmission measurement is performed in axial direction through the central volume of the tablet centre. However, in the subsequent analysis the measurement is considered to represent a function of nominal porosity for the entire volume of the tablet. This is only strictly valid in the absence of density distribution in the radial direction. Most tablets in the pharmaceutical industry are of biconvex geometry (Kadiri and Michrafy, 2013) and it is well-established that the punch geometry influences the pressure distribution in the tablet during compaction as well as during ejection (Kadiri and Michrafy, 2013;Mazel et al., 2015). Other than for flat-faced tablets this can result in density distribution for a range of common tablet shapes, such as in biconvex tablets (Sinka et al., 2004). The biconvex shape of the punch results in a pressure gradient in the radial direction from the centre to the edges. This may result in a lower pressure in the tablet centre compared to the edges during compaction and it has been shown that the punch curvature can lead to a density distribution (i.e. porosity distribution) in the radial direction of a biconvex tablet (Sinka et al., 2004;May et al., 2013).
THz-TDS has been suggested as a method suitable for application as a process analytical technology to measure and control tablet porosity (Bawuah and Peiponen, 2016;Bawuah et al., 2016c,b). In this study, we aim to provide a thorough analysis of the merits and limitations of this technique together with a discussion of the potential sources and magnitude of error and how they can be mitigated and controlled. In addition, a new concept based on the AB-EMA and the Wiener bounds effective medium approximation is presented to investigate changes in the pore shape with porosity as well as to evaluate the applicability of the AB-EMA on a set of test data.

The Anisotropic Bruggeman Effective Medium Approximation (AB-EMA)
The traditional Bruggeman model is a commonly used approach for EMA that was developed in the early 20th century by Bruggeman (1935). It describes the dielectric properties of a composite material (the so-called effective medium) by combining the individual dielectric properties, ε, of the volume fractions occupied by the isolated materials where the minor fraction, the inclusions, are assumed to be of spherical shape and the major fraction forms a continuous phase (Bruggeman, 1935;Bawuah et al., 2020). This formalism was later developed further by including a depolarisation factor into what is now known as the AB-EMA, where inclusions are spheroids with their polar axis aligned parallel to the propagation direction of the electromagnetic field (Jones and Friedman, 2000). Eq. 1 describes the AB-EMA model for a porous two component system composed of air and a solid material, denoted by the subscript s, where ε s and ε eff denote the complex dielectric permittivities of the solid fraction and the effective medium respectively and f is the porosity. Here it was assumed that the complex dielectric permittivity of air is equal to 1.
The depolarisation factor, L, describes the overall shape of the inclusions, which in case of porous pharmaceutical tablets are the tablet pores. It can be converted into the aspect ratio of the semi-axes of the spheroid using Eq. 2 where the polar axis, p, is parallel and the equatorial axes, q and r, are perpendicular to the wavevector of an electromagnetic wave, with q = r (Jones and Friedman, 2000). For p = q, L becomes 1/3 and the inclusion takes a spherical form (and Eq. 1 results in the traditional Bruggeman model). Applying Eq. 2 results in a one dimensional, needle shaped prolate for L = 1 and a two dimensional, flat oblate for L = 0.
For multiple solid components it was demonstrated that the dielectric properties of the different solid components can be combined into a single solid fraction with the average dielectric properties (Bawuah et al., 2020).
Most pharmaceutical materials exhibit relatively low absorption losses at terahertz frequencies. Therefore, it is often assumed that the complex part of the dielectric permittivity of the solid fraction can simply be neglected, resulting in ε s = n 2 s and ε eff = n eff 2 , where n s and n eff are the intrinsic refractive index of the combined solid fractions and the effective refractive index, respectively (Bawuah et al., 2020). Under these assumptions, Eq. 1 simplifies to Eq. 3 which is the common expression used in terahertz porosity analysis.

Bounds for L
The Wiener bounds are a common EMA developed by Wiener (1912). They describe the two extreme cases in which a composite material is formed either through serial or parallel alignment of the individual components, which are assumed to behave like planar plate capacitors. This results in an upper and lower bound for ε eff of the composite material, which again can be approximated in terms of n eff for weakly absorbing materials (Wiener, 1912;Markl et al., 2018a).
Eqns. 4 and 5 show the upper and lower Wiener bound for the general case of a multicomponent material with J solid components with the corresponding mass fractions x j and refractive indices n j .
In case of a two component system of air and one solid fraction, s, Eqns. 4 and 5 simplify to Eqns. 6 and 7.
Note that in this concept multiple solid fractions are again represented by one refractive index, n s , the same approach used in the previous section. For this case the Wiener bounds could in principle also be extended to their multicomponent equivalents (Eqns. 4 and 5). However, these do not conform with the AB-EMA model which only considers two phase compound materials and are therefore not further discussed in this article. We refer to an article by (Bawuah et al., 2016a) where these equations where utilised in the pharmaceutical setting using terahertz analysis. Fig. 1 illustrates the concept of the upper and lower Wiener bounds in THz-TDS measurements. The true value of n eff must fall between the extrema defined by the Wiener bounds for any given porosity (Bawuah et al., 2016a). The concept of the Wiener bounds is advantageous in cases where the microstructure of a material is unknown since no assumptions need to be made with regards to the pore shape (Bawuah et al., 2016a;Markl et al., 2018a,a). Further, the Wiener bounds can be applied independent of the sample shape.
Interestingly, the same expressions as in Eqns. 6 and 7 can be derived directly from the AB-EMA model. Assuming a depolarisation factor L = 0 (prolate, needle shaped pores) or 1 (flat, oblate pores) results in the upper and lower Wiener bound, respectively (see Fig. 1). It is therefore possible to convert the Wiener bounds into bounds of L in cases for which n s is unknown. For f = 0 the bounds converge to n s . The best estimate of n s in a porous material can therefore be achieved at the lowest porosity where the Wiener bounds are closest together. Eqns. 8 and 9 show these modified expressions of the Wiener bounds resulting in two margins for n s based on the sample with lowest porosity (f 1 ) and its effective refractive index (n eff, 1 ). It should be noted that the value of n s of a material is a physical constant, similar to the concept of a true density of a material.
Other than n eff in the Wiener bounds (Eqns. 6 and 7), the bounds in Eqns. 8 and 9 do not define values that n s can take dependent on the pore structure but rather describe the margins it must lie in between. The margins of n s can be interpreted as the refractive index that the solid fraction must have if the lowest porosity sample has pores (described by L 1 ) for which L 1 = 0 or 1 (or in the Wiener bounds formalism: serial or parallel aligned pores) to result in the measured value of n eff, 1 .
Replacing n s in Eqns. 3 with each margin results in upper and lower bounds (L u and L l ) for L: The expressions in Eqns. 8-11 can be used to estimate the values of n s and L when their true values are unknown.
An immediate release formulation of ibuprofen (BLD Pharmatech, Shanghai, China) was made with dose strength of 10% w/w. The Fig. 1. Illustration of the Wiener bounds following the conventional capacitor model and the AB-EMA inclusion shape. The spheroid represents the orientation of the air pore in the dielectric medium of the solid phase relative to the propagation direction of the electromagnetic field of the transmitted terahertz pulse. formulation was composed of common excipients, MCC, lactose anhydrous (Supertab21AN, DFE pharma, Goch, Germany), croscarmellose sodium (CCS, DuPont Nutrition, Wilmington DE, USA) and magnesium stearate (Fisher Scientific, Fair Lawn NJ, USA).

Tableting
The powders were directly compressed with 10 mm round flat-faced and biconvex tooling on a Huxley Bertram 50 compaction simulator (Huxley Bertram Engineering, Cambridge, UK) using a symmetrical compaction profile. The punch depth and radius of curvature of the curved punches for biconvex tablets was 1.05 mm and 12 mm, respectively. The compaction profile was based on a Fette 102i with a rotation speed of 20 rpm corresponding to a die speed of 293mm s − 1 . Each main compaction was proceeded by a pre-compression step with a compression gap of 7.5 mm. The pressure during the main compression was in the range of 5160 MPa. Tablet mass was kept constant at 300 mg, for flat-faced, and 400 mg for biconvex tablets. The material for each tablet was individually weighed, manually filled into the die and compressed.

Porosity determination
The tablet porosity of flat-faced and biconvex tablets was determined based on Eqns. 12 and 13, respectively.
where M is the tablet mass, D is the tablet diameter, H is the tablet height, h is the punch depth of curved punches, C is the radius of curvature of curved punches, and ρ true is the true density with values of 1.527 kg m − 3 for the MCC blend and the lactose powder, 1.439 kg m − 3 for the ibuprofen blend and 1.490 kg m − 3 for the starch powder.

THz-TDS
THz-TDS measurements were acquired in transmission geometry with the beam propagating through the entire tablet using a Terapulse 4000 spectrometer (Teraview Ltd., Cambridge, UK). The measurement cell was purged with dry nitrogen gas (relative humidity <1%) to avoid any pertubation by water vapour. For each measurement 20 waveforms were averaged. An empty measurement cell was taken as reference and the reference signal was acquired before each tablet measurement.

Data analysis
Data analysis and simulations were carried out using MATLAB 2020a (MathWorks, Natick, MA, USA). Extraction of n eff was performed by taking the ratio of the Fourier transform of the time-domain data for each sample and reference pulse. A more detailed description of the extraction from the optical constants from the time-domain can be found in Jepsen and Fischer (2005). For the data analysis n eff was averaged in the frequency range of 0.4 -0.8 THz. The AB-EMA fit was carried out by applying Eq. 3 to the data by varying L between 0 and 1 in increments of 0.001. The fit with the smallest standard deviation in n s at different porosities was considered best.

Data simulation
Two phase powder compacts of air and solid fraction were simulated based on the AB-EMA (Eq. 3).

Effect of material absorption
To simulate the effect of material absorption on the AB-EMA model samples were simulated with α eff ranging from 0 cm − 1 to 50 cm − 1 in increments of 1 cm − 1 , n eff ranging from 1 to 5 in increments of 0.1 and porosities from 0.05 to 0.25 in increments of 0.05. For each sample κ eff was calculated by using the expression κ eff = ( cα eff )/ (4πν) where c is the speed of light. The frequency ν was set to 0.2 and 1 THz. Eq. 1 was then used to calculate n s with and without accounting for absorption. For the latter case Eq. 1 results in Eq. 3. The pore shape was evaluated based on Eq. 2 with a p : q ratio of 5, 1 and 1/5 resulting in L values of ≈0.053, ≈ 0.33 and ≈0.75.

Variations in the pore shape
To simulate a change in L with porosity, a linear model L(f) = a 1 f + a 2 was assumed. The simulated porosity ranged from 0.05 to 0.5 in increments of 0.05. The value for n s was based on the current in-house estimate of n s of the common excipient MCC which is at 1.86. Different sets of tablets were simulated by stepwise changing L from 0 to 1 and 1 to 0 for the highest and lowest porosity, respectively in increments of 0.02. Porosities in between were adjusted accordingly. This resulted in the gradient, a 1 , taking on values from -2.22 to 2.22 in increments of 0.0889 and resulted in the averaged L over all porosities being constant at 0.5.

Density distributions of biconvex tablets
To simulate biconvex tablets with different density distributions, the nominal porosity of the tablets ranged from 0.05 to 0.5 in increments of 0.05. The value for n s was 1.86 as described above. The centre porosity was set 10% higher and lower than the nominal porosity to simulate a lower and higher density in the tablet centre. The corresponding n eff values were then calculated for the different centre porosities based on Eq. 3 with L ≈ 0.33 (spherical pores).

Limitations of the AB-EMA
For most materials n s cannot be determined directly since a zeroporosity sample is experimentally inaccessible. Instead, n s is determined by an extrapolation from the regression analysis of n eff over the experimentally accessible porosity range using the AB-EMA (Bawuah et al., 2020).
The correct extraction of n s and L is of critical importance since they define the AB-EMA fitting model, and therefore the correlation function, used for porosity prediction. Furthermore, knowledge of n s is required for more advanced pore shape analyses based on the Wiener bound model (Bawuah et al., 2020). In this study the estimated values of the AB-EMA fit for n s and L will be referred to as n s, fit and L fit . The range of α eff was based on a non-absorbing material with α eff = 0 and a "worst case" scenario in which α eff was based on a thin tablet of 1 mm thickness with a transmittance of 1%. The porosity range of 0.05 to 0.25% was based on the typical values used for pharmaceutical tablets. For n eff , common literature values were considered for setting the range (Lu et al., 2020). The two frequencies investigated are typical limits of the dynamic range of the spectrometer used in our lab when measuring entire tablets in transmission. For L Eq. 2 was utilised at a p : q ratio of 5 (prolate), 1 (spherical) and 1/5 (oblate).

Effect of material absorption on the AB-EMA
The data shown in Fig. 2 demonstrate that the relative error is highly dependent on n eff , α eff , L, the porosity and the frequency. However, it should be noted that in any material absorption decreases with increasing porosities. The magnitude of the relative error was found to cover the range of 1 × 10 − 0.7 . A general trend towards larger errors at higher absorption and lower values for n eff can be observed. However, in some simulated sets intermediate values of n eff lead to "islands" of higher errors (e.g. at L ≈ 0.33, f = 0.25, ν = 0.2 THz). An increase in the effect of absorption can be observed at lower frequencies which is due to the frequency dependence of κ eff (see Section 2). This should be kept in mind during data analysis since clearly a critical threshold is reached earlier if optical constants are extracted at lower frequencies. The relative error additionally shows a remarkable dependence on the pore shape, described by L. Interestingly, under the investigated pore shapes, spherical pores (L ≈ 0.33) resulted in the smallest error. For nonspherical a more oblate shape (L ≈ 0.75) seems to result in higher errors compared to more prolate pores (L ≈ 0.053). Data extraction of n eff typically results in errors on the order of 1 × 10 − 4 to 1 × 10 − 3 (see Supplements Tables 1 to 6). Therefore, in most cases the error caused by the simplified model will be less than that caused by the measurement system itself and the effect due to absorption will be negligible. For comparison, tablets of the typical pharmaceutical excipients MCC, lactose and starch as well as an ibuprofen formulation are shown for which small errors of magnitudes <10 − 4 can be expected. Please note that the pore shape of these tablets was not investigated and therefore the materials are displayed for all three investigated L values.
In non-crystalline materials, such as the majority of pharmaceutical excipients, the absorption spectrum exhibits monotonously increasing absorption without any sharp spectral features. In fact, this increase in absorption constitutes the rising flank of the peak due to the vibrational density of states (Sibik et al., 2014). This peak is wider than the spectral bandwidth of any typical THz-TDS instrument that is currently commercially available. Such spectrometers commonly access a range of about 0.1 THz to 3 THz, and in some cases up to about 7 THz, but the dynamic range of all such instruments is decreasing rapidly with frequency and hence for measurements of any tablet sample the usable bandwidth is limited to a maximum of about 2 THz. Therefore, one might expect higher errors at higher frequencies due to the stronger absorption as well as the reduced dynamic range. However, due to the frequency dependency of κ eff the effect of increasing absorption is less pronounced. This can also be observed in Fig. 2: although absorption is significantly higher at 1 THz the expected error remains at lower values.
In conclusion, neglecting material absorption during data analysis and application of the AB-EMA is justified in most cases, and in particular for the materials investigated here. The choice of materials used in this study was aimed to represent a wide range of material and structure properties that are typically encountered in pharmaceutical solid dosage forms. In general it is recommended to chose a suitable frequency window for the AB-EMA analysis where the tablet matrix exhibits relatively low absorption. This ensures the validity of the simplified AB-EMA model and further ensures a constant refractive index over a large frequency range in cases where high material absorption is due to a resonance peak in the spectrum. However, extra care must be taken to ensure the assumption is justified when applying the method to materials known to form highly anisotropic pore structures upon compaction as well as with materials with high absorption at low frequencies, as some materials may exhibit values for α eff above the considered range. In cases where the validity of these assumptions is uncertain we recommend to extend the AB-EMA model by including the absorption coefficient and using the complex refractive index instead of its real part only. This will effectively eliminate the errors associated with absorption and the AB-EMA model. Since THz-TDS allows for the direct extraction of both parts of the complex refractive index without requiring any additional tools, such as chemometrics or advanced mathematical models such as the Kramers Kronig analysis, this will not affect the experimental setup but only the complexity of the applied model. However, for tablets comprised of such strongly absorbing materials the deterioration of the signal-to-noise-ratio for a measurement in transmission is likely to be a more practical limitation that may prevent application of this method altogether.

Effect of variations in the pore shape on the AB-EMA
To investigate how porosity dependent changes in the pore shape, described by L in Eq. 3, affect the accuracy of the AB-EMA fit, powder compacts were simulated as described in Section 3.6. It should be noted that for gradients above 2.04 the algorithm was unable to fit the data. For a gradient of 2.04 the corresponding spheroids have a p : q ratio of 6.0 and 0.026, at the highest and lowest porosity, respectively (Eq. 2). Such drastic changes in pore shape are unlikely to occur in practice and therefore this was not investigated further in this study. Fig. 3 illustrates the results for an example set with the simulation parameters, a 1 = 0.5 and a 2 = 0.275 (Set 1). Depicted are the simulated values of n eff , the different AB-EMA models for each porosity used in the simulation and the AB-EMA fitting function. Just by visual inspection it is clear that the AB-EMA fit was unable to follow the simulated data points and further did not result in a fitting function similar to the simulated AB-EMA models. Fig. 4 illustrates how the adjusted R 2 and the root-mean squared error (RMSE) for the porosity are affected by the gradient (simulation parameter a 1 ) of L. The RMSE was calculated for the independent variable (the porosity) instead of the dependent variable (n eff ). This was considered more practical since for tablet characterisation n eff is used to estimate the porosity. For the set with a 1 = 0 both parameters reached a minimum close to zero. This was to be expected since L remained constant for this set and could therefore be described by Eq. 3. With higher absolute gradients the fit was impaired gradually. However, the effect on the goodness of fit based on the adjusted R 2 was marginal and even for strong variations in the pore shape the fit was able to describe the data with values above 0.99 at all gradients.
Based on the analysis of the RMSE, deviations in the porosity prediction can be expected to fluctuate by a maximum of around 5% for the theoretically most drastic changes in pore shape. For realistic changes that are more likely to be encountered in tablet processing much lower deviations are to be expected. However, the error was not consistent throughout different porosities as illustrated by the residual surface plot (Fig. 4). Larger residuals were found for the centre porosity and the high and low porosity ends. Therefore, larger deviations between the true and predicted porosity are to be expected at these porosities.
The threshold up to which error a method for porosity prediction is useful will very much depend on the final formulation and the demands on quality and performance of the product. To the best of the authors' knowledge, guidelines currently do not provide specifications regarding tablet porosity. Since extreme gradients in L with porosity are unlikely to occur in practice, large deviations in real sets of tablets due to changes in the pore shape are not considered a practical problem. Due to the porosity dependent distribution of the residuals, one should consider an arbitrary fitting model instead of the AB-EMA fit when large deviations in the pore shape can be expected.
The fitting method for the AB-EMA model selects the value of L fit for which applying Eq. 3 results in the least variation of n s, fit at the different porosities. The final value of n s, fit is then calculated by averaging. Fig. 5 shows a plot of the relative error of n s, fit at all porosities, as well as its average, as a function of the gradient a 1 . For sets where L varied, the gradient strongly affected the accuracy of the fit. For positive gradients an overestimation of n s, fit was observed whereas the opposite was the case for negative gradients.
For relatively modest absolute gradients (m > | 0.6|) relative errors in excess of 2.5%, were determined even without considering any additional experimental errors. In this context, it is interesting to note that the estimates of n s for the upper and lower porosity limit were least affected and their relative error remained close to 0. For L fit no single reference value existed since a range of values were used in the simulation. L fit is considered to be a rough approximation for the description of the pore shape and thus the range of L values that were chosen in the simulation (Markl et al., 2018b).
The bottom panel of Fig. 5 shows L fit together with the maximum, L max , and minimum, L min , L values of the simulated data. As in the top panel of Fig. 5, the AB-EMA fit resulted in a perfect estimation of L for the zero-gradient set. For most other gradients L fit lay outside the range of L values used in the simulation (represented by the grey area in Fig. 5). This result is quite unexpected and disproves the aforementioned hypothesis of L fit as a rough approximation of L. Instead L fit did not describe the pore shape of most of the simulated tablet sets since it lay outside the range of L values used in most cases. As for the relative error of L fit , negative and positive gradients resulted in an under-and overestimation, respectively. L fit was always closest to the highest porosity sample for which at negative gradients L min and for positive gradients L max is the corresponding L value.
Given these results, it is important to be mindful of the limitations of this approach when adapting the AB-EMA method to data sets where the pore shape, and therefore L, cannot be assumed as constant. The tableting process is anisotropic in its nature and deviations in the anisotropy with the pressure can often be expected (Porion et al., 2010). The model of using a linear change in L with porosity that was used for this study is only a simplified model for illustrative purposes. In reality more complicated, and indeed formulation and process dependent, deformation behaviour can be assumed. Despite the over simplification of our model the data illustrate that false estimates of n s and L can be expected even for small changes (i.e. small gradients) in the pore shape upon compaction to different porosities.
Based on the data shown in Fig. 5 one might argue to base the estimate of L on the sample with the highest porosity since for this sample the error of n s, fit and L fit is the smallest. However, this section demonstrated that the AB-EMA model is inadequate to fully describe the pore shape based on THz-TDS measurements and and hence this approach is not recommended.
For porosity prediction the effect of variations in L with the pore shape is not considered a critical issue given that the adjusted R 2 values are close to 1 and the RMSE values are very low. Therefore, the AB-EMA is still considered the most suitable approach for porosity prediction. In cases where the AB-EMA fit of a calibration set results in large residuals one can consider the application of other fitting functions which better describe the formulations changes of n eff with porosity in order to get a more rigorous porosity estimation. However, a physically meaningful model, such as the ones based on the formalism of the effective medium approximation, will enable us to learn about the formulation mechanistically beyond simple porosity prediction and is therefore more desirable over applying non-physical fitting functions. Nonetheless, our results indicate that a good correlation of the model alone must not be assumed to indicate the correct extraction of n s and L. A more rigorous method to estimate the applicability of the AB-EMA to a set of data will be introduced in Section 4.2. Fig. 6 illustrates the THz-TDS setup for biconvex tablets. The terahertz beam propagates through the centre of the tablet from the top to the bottom tablet face. It therefore only measures n eff of the centre porosity. In contrast, the nominal porosity of the entire tablet is calculated using the tablet's dimensions and mass and therefore represents the averaged porosity of the sample assuming even distribution of density. In practice the porosity of a sample at its centre is unknown.

Effect of density distribution in biconvex tablets on the pore shape prediction
To investigate the effect of density distribution on the AB-EMA, tablets were simulated as described in Section 3.6 with n s = 1.86 and L ≈ 0.33. Even density distribution and a 10% higher and lower centre porosity compared to the nominal porosity were simulated in three different sets. Density distribution was loosely based on the results of the finite element analysis method for biconvex MCC tablets as reported by Sinka et al. (2004), where differences up to 30% between the tablet centre and the nominal porosity were reported. Based on these results

Fig. 5.
Error analysis of the AB-EMA fitting parameters as a function of the gradient. The relative error of n s, fit as a function of the gradient as a mean for all porosities (black, dot) as well as for each individual porosity shown on the top. L fit (black, dot), L max (green, dot), and L min (red, dot) as well as all L values in between used in the simulation (grey, dashed) as a function of the gradient (bottom). The grey area illustrates the range of L used for each gradient. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 6. Schematic of the geometry of the beam propagation in a terahertz measurement setup illustrated for a biconvex tablet. The terahertz beam penetrates the centre of the tablet from one face to the other and thereby only probes the centre porosity. In typical setups the beam waist will have dimensions of a couple of millimetres. However, in biconvex tablets, density distribution (illustrated by different shades of blue) might differ from the nominal porosity estimated by tablet weight and dimensions. Please note that this illustration is only accurate for single wavelength beams as the beam width of a Gaussian beam will be frequency dependent. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the simulated deviation of 10% was considered realistic. Fig. 7 shows a plot of n eff and the AB-EMA fits as a function of nominal porosity. Compared to t he set with even density distribution, n eff decreased and increased with higher and lower centre porosity, respectively. This was a consequence of the difference between the nominal porosity and the porosity extracted from the measurement of the tablet centre, which resulted in a corresponding change of n eff . Naturally, the AB-EMA fit of the evenly distributed sample led to a perfect correlation with the adjusted R 2 being equal to 1. For samples with higher and lower centre porosity the correlation was only marginally impaired with adjusted R 2 values of 0.9993 and 0.9997. L fit was found at 0.498 and 0.160 for the sets with higher and lower centre porosity, respectively. Thus, density distribution resulted in a false estimation of L, which was fixed to ≈0.33 in all samples. By using Eq. 2, this correspond to a p : q ratio of approximately 1 : 2 and 2 : 1 for the high and low centre porosity, respectively. Without knowledge of the true cause both could easily be mistaken for plausible pore shapes.
With new insights gained from the analysis of simulated data, nonsimulated, flat-faced and biconvex MCC tablets were investigated. The tablet properties are listed in Table 5 and Table 6 in the Supplementary Materials. Fig. 8 shows a plot of n eff as a function of porosity for the two types of tablets. Averaging in the frequency range of 0.4 -0.8 THz was performed to account for any small fluctuations in n eff as a function of frequency. The frequency range was based on a previously suggested range for porosity measurements of MCC tablets (Bawuah et al., 2020). For the biconvex tablets lower values of n eff were observed compared to the flat-faced samples. This suggests a higher centre porosity compared to the nominal porosity. For the AB-EMA fit a good correlation of n eff and porosity was found with slightly better results for the flat-faced tablets. Adjusted R 2 values were at 0.9993 and 0.9874 for flat-faced and biconvex tablets, respectively. For the biconvex tablets L fit reached a higher value of 0.415 compared to 0.357 for flat-faced tablets.
Based on the previously gained results for the simulated data it is suggested that the higher value of L fit in the biconvex tablets originated from density distribution in this set. In this case an assumption of equal centre and nominal porosity would result in an apparent change in the pore shape compared to flat-faced tablets when using the AB-EMA. Based on the simulation, the higher value of L fit corresponds to a higher centre porosity. This is well in line with the lower values found for n eff in this set compared to the flat-faced tablets.
Several studies suggest that biconvex tablets often show higher porosity in the centre compared to the edges (Sinka et al., 2004;May et al., 2013;Diarra et al., 2015). In contrast, flat-faced tablets are not expected to show significant density distribution in the radial direction, if any . Therefore, the set of flat-faced tablets can be used for the accurate extraction of the material property, n s , and thus to reliably predict the centre porosity of biconvex tablets.
It should be noted that for both, biconvex and flat-faced tablets, density distribution in the axial direction, that is the direction the force is applied in during compaction, may very well be present. The porosity obtained by THz-TDS is always an average measurement that is not able to resolve any axial density distribution. It is well-known that the vast majority of axial density anisotropy is removed during elastic recovery in the unloading step within the die of the tablet press and the effects due to any residual remaining strain in the tablet is the subject of current research.
Based on the terahertz analysis, an average ratio between the porosity measured at the tablet centre and the nominal porosity of 0.875 was found for the biconvex tablets investigated here. As already concluded from the lower n eff values for biconvex tablets, the centre porosity was therefore higher compared to the nominal porosity. Considering the similar density distribution in the real tablets and the simulated tablets with lower centre porosity one would expect a higher deviation in L fit between the biconvex and flat-faced tablets. This could be rooted in the high variation of density distribution between tablets of different porosities (See Supplements Fig.12). Taking this into account the experimental and simulated deviations in L fit are considered to be in good agreement by the authors.
Data analysis of simulated data demonstrated that density distribution in biconvex tablets highly affects the AB-EMA methodology when the nominal porosity of the tablet is considered equal to the centre porosity. In this case the AB-EMA fit results in an inaccurate estimation of L fit . Similar as the variation in pore shape discussed in Section 4.1.2, the goodness of fit was not impaired significantly by density distribution. This again creates a situation where one might not be able to establish whether the extracted parameters are accurate based on the Fig. 7. n eff and AB-EMA fit as a function of nominal porosity for simulated tablets with density distribution. The results are depicted for tablets with even density distribution (blue) and a centre with 10% higher (red) and lower (green) than average porosity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)  Table 5 and Table 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) ability of the AB-EMA to correlate n eff and porosity. Previously, a change in pore shape with tablet geometry was reported by applying the AB-EMA on THz-TDS measurements (Bawuah et al., 2019). Based on the results discussed in this section it is suggested that this may be rooted in the density distribution of biconvex tablets rather than an actual change in pore shape. This is consistent with the density distribution suggested for the here reported biconvex MCC tablets.
In conclusion, it is recommended to perform all calibration measurements using flat-faced tablets rather than biconvex or other shapes that may result in an anisotropic density distribution. The results obtained from the flat-faced tablets can then be directly used together with the AB-EMA method to accurately measure the centre porosity of biconvex tablets, a metric which cannot be determined easily by other methods. The tablet centre porosity may in many cases constitute a key critical quality attribute of an immediate release tablet with regards to disintegration due to its commonly higher porosity compared to the rest of the tablet matrix.

Extensions of the AB-EMA model to describe variations in the pore shape
Section 4.1.2 demonstrated how tablet anisotropy affects the AB-EMA and the extraction of n s and L. Although the AB-EMA resulted in a strong correlation, it was unable to extract the true physical properties of the material when pore shape variations or density distribution were simulated. Therefore, the accuracy of the extracted parameters cannot be assessed based on the goodness of fit alone. This section discusses alternative models to evaluate the applicability of the AB-EMA fit and to investigate L as a function of porosity.

Simulated data
To evaluate the new approach based on the concept of Wiener bounds, three sets of simulated data were used with corresponding linear models of a 1 = 0.5 and a 2 = 0.275 (Set 1), a 1 = − 0.5 and a 2 = 0.55 (Set 2), and a 1 = 0 and a 2 = 0.5 (Set 3). Fig. 9 depicts L l and L u for the three sets. The bounds were chosen based on the sample with lowest porosity and thus information got lost for this sample and the upper and lower bound resulted in values of 0 and 1 by default.
With increasing porosity, the bounds approach each other and their tangents ultimately converge parallel to one another and the simulated gradient of L. For Set 3 this results in a convergence towards a constant value of L, whereas for Sets 1 and 2 a close to linear in-and decrease, respectively, towards the high porosity end can be observed. Fig. 9 also depicts the values of L fit of the AB-EMA fit for the three sets. In Section 4.1.2 it was discussed that variations in the pore shape can result in a false estimate of L fit . This can also be seen when comparing L fit with the bound progression in Sets 1 and 2. There is a clear deviation between the constant progression in L fit and the gradient for the bounds at higher porosities. Furthermore, for Set 2 L fit remained below the bounds for porosities between 0.1 and 0.45. This is not possible since the bounds display the two extremes in which L must lie in.
Therefore, the bounds can be used as a simple method to validate the AB-EMA fit. If the bounds converge towards a steep gradient, or if L fit lies outside the calculated bounds at certain porosities, the AB-EMA fit is invalid and other methods need to be applied. No prior knowledge about the deformation behaviour of the material needs to be obtained. The THz-TDS data can be directly used to estimate the validity of the AB-EMA fit to avoid errors associated with a change in the pore shape (as discussed in Section 4.1.2) and to chose alternative models for a better description of the investigated tablet set, one of which will be introduced in Section 4.2.2. Furthermore, for high porosities a good approximation of the pore shape can be made due to the proximity of the bounds. Lastly, bounds of multiple sets with different pore shapes can be used to limit the range of possible L values.

Comparison of different size fractions of MCC applying the bounds for L
Skelbaek -Pedersen et al. (2020a) previously investigated the effect of fragmentation during tableting on the water ingress kinetics using terahertz pulsed imaging. Their data was reused in this section to evaluate the newly introduced concept of upper and lower bounds for L and the effect of deformation on the AB-EMA on non-simulated tablets. Two sets of tablets compacted with different MCC particle size fractions in the range of <125 μm and 355 μm to 500 μm were analysed by THz-TDS.
The materials were compressed into round, flat-faced tablets at different pressures resulting in porosities of 0.05 to 0.3. The powder mass was adjusted accordingly to retain a constant tablet thickness. For further discussion, the two powders will be referred to as small and large size fractions for particle sizes of <125 μm and 355 μm to 500 μm, respectively. Fig. 10 depicts n eff and the corresponding AB-EMA fit as a function of porosity for tablets of both size fractions. For the larger size fraction higher values of n eff were measured. Both size fractions are composed of the same material and therefore, n s must not change between size fractions. The fitting method was adjusted to result in the best model for which n s, fit takes on the same value for both size fractions. The AB-EMA fit estimated n s, fit = 1.867 and L fit at 0.259 and 0.339 respectively for the large and small size fraction.
This means that the AB-EMA predicted more oblate shaped pores for tablets of the small compared to the large size fraction. To evaluate the accuracy of the fit, L l , L u , and L fit are plotted as a function of porosity in the bottom panel of Fig. 10. Measurements of tablets samples prepared from the small size fraction powder resulted in higher values of L l and L u . Thus, similar as for L fit the bounds suggest a more oblate pore shape for the small size fraction. The values of L fit fell outside the bounds for porosities between 0.20 -0.30 for the large size fraction and porosities at 0.10; 0.20; 0.25; 0.33 for the small size fraction.
For these samples the extracted values of L fit must be erroneous and the AB-EMA model thus failed to describe the data. The progression of the bounds at high porosities suggests a negative gradient of L for both tablet sets. Such a change in L would likely result from an incorrect estimate of L by the AB-EMA fit as discussed in Section 4.1.2. Furthermore, based on Fig. 4 the negative gradient should result in an underestimate of L fit compared to the range of true L values of the pores. This is in accordance with values of L fit falling below the bounds at some porosities for the two types of tablets. It can therefore be concluded that the values of L fit must indeed be underestimated at these porosities. The results suggest that the AB-EMA method is inadequate to describe the two sets of tablets because of pore shape changes in these tablets under the conditions studied.
To adjust the AB-EMA and allow for changes in the pore structure with the porosity to be accommodated, the L term can be replaced by a porosity dependent expression. It was found that replacing L in Eq. 3 with a simple linear model of the form L(f) = a 1 f + a 2 can describe the data well. It is helpful to highlight that replacing L with this model still allows for a constant value of L in cases where the gradient, a 1 = 0. The modified AB-EMA (mAB-EMA) model is shown in Fig. 10 as a function of porosity. The mAB-EMA model estimated higher values of n s at 1.869. It resulted in slightly higher adjusted R 2 values of 0.9942 and 0.9991 for the large and small size fraction compared to the AB-EMA (adjusted R 2 at 0.9936 and 0.9984). L fit of the mAB-EMA lay inside the bounds over the entire range of porosities (Fig. 10). The model resulted in a clear gradient of − 0.141 and -− 0.289 of L fit as a function of porosity for the large and small size fraction, respectively. This further supports a change in pore shape towards more oblate pores at lower porosities as already suggested by the progression of the bounds. Additionally, it resulted in larger values of L fit for the smaller size fraction and therefore suggests more oblate pores compared to the large size fraction, in line with the other models.
Without any further data processing the data shown in Fig. 10 already suggest a difference in tablet microstructure in the two size fractions based on the different values of n eff at a specific porosity. For tablets of the same material with the same microstructure n eff should result in the same value. However, while the work by SkelbÃ¦k-Pedersen et al. did report a difference in microstructure between tablets of different MCC size fractions (Skelbaek-Pedersen et al., 2020a, 2020b, 2020c, in their study on measuring water ingress into the porous tablet matrix it was not explored specifically that THz-TDS is able to measure this difference in structure but the work focused on the liquid transport kinetics (Skelbaek-Pedersen et al., 2020a). Further analysis of the data show that the AB-EMA model, L l and L u , and the mAB-EMA both suggest more oblate shaped pores, based on higher L values, for the tablets prepared from the small compared to the large size fraction.
It was furthermore reported that an increase in particle size results in a higher fragmentation degree during tableting (Skelbaek-Pedersen et al., 2020c). Small size fractions are therefore expected to show less fragmentation. This was demonstrated on the same size fractions of MCC as used in this study. If less of the energy during compaction is converted for two sets of MCC tablets with initial particle sizes of <125 μm (blue) and 355 μm to 500 μm (red). Error bars indicate the standard deviation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) into particle fragmentation, plastic deformation must occur to a higher degree for tablets of the same final porosity given the amount of elastic deformation is materials dependent and therefore constant. For MCC it is well-known that plastic deformation is the major deformation route (Skelbaek-Pedersen et al., 2020c) and therefore tableting might result in a higher degree of plastic deformation in smaller size fractions due to a lesser degree of fragmentation.
Such behaviour would be in line with the pore shape analysis presented in this section. The higher degree of plastic deformation would result in relatively more flatter particles in the resulting tablets. Since pore shape is considered the indirect result of the shape of the surrounding particles this would in turn result in a flattening of pore shape as well. Wu et al. (2008) reported the same for cubic starch compacts based on image analysis methodology. Upon compression particle deformation resulted in "pancake" like structures. Pores were formed by stacking of the particles and in turn exhibited high anisotropy (Wu et al., 2008).
In the AB-EMA model this would correspond to more oblate pores with L values closer to 1 as observed in this study. An effect of plastic deformation on L can be further suggested based on the negative gradient of L fit obtained by the mAB-EMA as well as indicated by L l and L u . This gradient is likely to be the result of particle deformation upon higher compression pressures. Therefore, both models, the bounds for L and the mAB-EMA, predicted a change in the pore shape well in line with theoretical considerations regarding plastic deformation of the two tablet types.

Effect of density distribution on the bounds for L
In Section 4.1.3 the effect of density distribution in biconvex tablets on the AB-EMA method was discussed. The same sets of simulated and real biconvex and flat-faced tablets are re-analysed for this section. To evaluate the effect of density distribution on the bounds, L l and L u , of the simulated tablets were plotted as a function of porosity in Fig. 11.
For the three sets the bounds progressed differently. This suggests a difference in the porosity dependent pore shape. However, L was kept constant at ≈0.33 in the simulation. The difference in the observed trend is therefore purely the result of density distribution and the false assumption of equal centre and nominal porosity. Here higher centre porosities (as expected for most biconvex tablets, (Sinka et al., 2004;May et al., 2013;Diarra et al., 2015)) resulted in a strong apparent dependency towards more oblate pores at lower porosities. Fig. 11 illustrates the same analysis on non-simulated biconvex and flat-faced MCC tablets. The porosity range of the two tablet sets, while being similar, was not the same. This resulted in the bounds for the flatfaced tablets being further apart since in this set the lower end of the porosity range was higher compared to the biconvex tablets. For the flatfaced tablets, no clear trend in L l and L u as a function of porosity was observed. For the biconvex tablets, L l and L u indicate a strong trend towards higher values at lower porosities. Based on this, one would expect the AB-EMA to be unable to represent the data, which is consistent with the results in Section 4.1.3.
Since for biconvex tablets density distribution is common, one should keep this in mind when attributing the trend in L l and L u to a change in pore shape. Based on the results for the simulated data it is suggested that the strong dependence is in fact the result of the higher centre porosity estimated in Section 4.1.3. This itself is not necessarily a disadvantage when L l and L u are used to test the applicability of the AB-EMA fit. Steep gradients of the bounds at higher porosities indicate that the AB-EMA cannot be applied for the analysed set. It does not matter whether this is due to density distribution or pore shape changes. Furthermore, it can be useful when investigating density distribution in biconvex tablets with THz-TDS. To differentiate between the effect of density distribution and pore shape on L l and L u the results between a flat-faced calibration set and the biconvex analysis set can be used. As demonstrated in this section a difference in n eff of flat-faced and biconvex tablets at the same nominal porosity indicates density distribution in the radial direction. When the centre porosity is expected to differ from the nominal porosity, L l and L u should not be used to investigate the pore shape of biconvex tablets.

Conclusion
Different error sources of the AB-EMA for the extraction of porosity and pore shape of pharmaceutical tablets with THz-TDS were investigated. Absorption was found to have a negligible impact on the accuracy of the model. Pore shape changes with porosity and density distribution did both not significantly impair the quality of the AB-EMA fit for porosity estimation of pharmaceutical tablets. However, density distribution and pore shape changes highly influenced the accuracy of the estimation of n s and L. Further, the only marginal impairment of the goodness of fit increases the danger of misinterpreting the extracted parameters.
An alternative model was suggested to convert the Wiener bounds into bounds of L. It was demonstrated that the model can be used to Fig. 11. Evaluation of L l and L u of biconvex tablets. The bounds are depicted for simulated biconvex tablets without density distribution (blue) and a higher (red) and lower (green) centre porosity (left) as well as for real flat-faced (blue) and biconvex (red) MCC tablets (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) evaluate i) whether changes in the pore shape with porosity occur and ii) the applicability of the AB-EMA fit. Furthermore, a modified AB-EMA (mAB-EMA) model was introduced allowing for a linear change in L with porosity. The mAB-EMA model was found useful as a correlation function as well as for pore shape analysis.
Lastly, density distribution was found to affect the newly introduced bound model for L. However, it was suggested that for such cases the model is still useful for evaluating density distribution although the physical description of pore shape is lost.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.