An Efficient Signal to Noise Approximation for Eccentric Inspiraling Binaries

Eccentricity has emerged as a potentially useful tool for helping to identify the origin of black hole mergers. However, owing to the large number of harmonics required to compute the amplitude of an eccentric signal, eccentric templates can be computationally very expensive, making statistical analyses to distinguish distributions from different formation channels very challenging. In this paper, we outline a method for estimating the signal-to-noise ratio for inspiraling binaries at lower frequencies such as those proposed for LISA and DECIGO. Our approximation can be useful more generally for any quasi-periodic sources. We argue that surprisingly, the signal-to-noise ratio evaluated at or near the peak frequency (of the power) is well approximated by using a constant noise curve, even if in reality the noise strain has power law dependence. We furthermore improve this initial estimate over our previous calculation to allow for frequency-dependence in the noise to expand the range of eccentricity and frequency over which our approximation applies. We show how to apply this method to get an answer accurate to within a factor of two over almost the entire projected observable frequency range. We emphasize this method is not a replacement for detailed signal processing. The utility lies chiefly in identifying theoretically useful discriminators among different populations and providing fairly accurate estimates for how well they should work. This approximation can furthermore be useful for narrowing down parameter ranges in a computationally economical way when events are observed. We furthermore show a distinctive way to identify events with extremely high eccentricity where the signal is enhanced relative to na\"ive expectations on the high frequency end.

However, the gravitational waves (GWs) radiated from inspiraling eccentric binaries possess multiple harmonics and the width of the spectrum increases with the eccentricity [22]. This complicates the calculation of the signal-to-noise ratio (SNR) as one has to sum over those many harmonics. The spread of frequencies makes the calculation of the average SNR more challenging in that one can no longer just calculate the signal and divide it by a fixed noise. This is especially problematic for large eccentricity since one needs to evaluate many Bessel functions at large orders, which is computationally expensive. Therefore it would be good to have a formula relating the SNR of an eccentric binary to the SNR of a more tractable system.
In this paper we elaborate and extend the simplified calculation for the SNR of eccentric binaries of [23], where two of us showed how to relate the SNR of an eccentric binary to that of a circular binary with the same masses and peak frequency f p , the frequency of the peak harmonic in the GW power spectrum. We proposed an approximating formula that assumed a constant noise curve, using both an analytical argument at large e and a numerical fit. Treating the SNR (f p , e) as a function of the peak frequency f p and the eccentricity e, we showed (f p , e)/ (f p , e = 0) ∝ (1 − e 2 ) 3/4 (1) for large e. We then provided a full formula working for all e ∈ [0, 1) through a numerical fit. For noise curves that decrease with frequency (or ones that decrease over the range of the signal), approximating the noise curve as a constant yields a partial cancellation between the frequencies that are lower and higher than the peak frequency. The corrections to our approximation would appear in the second derivative, which means it depends on the spread of frequencies divided by the peak frequency. So long as eccentricity is sufficiently far from unity, this will be a small(ish) number. Near the turnaround of the noise curve, there is less variation in noise from the get-go so the approximation there also works very well. Here we show that for eccentricities less than 0.9, this approximation works rather well (at least at the level of 50% accuracy and often better). This accuracy should be adequate for distinguishing among population distributions.
However, the partial cancellation does not necessarily work for the rising part of the noise curve because of the 1/n factor in the amplitude of the n'th harmonic, which manifests itself as a secondary peak in the spectrum at low frequency where the noise can be relatively low, as can be seen in Fig. 2. We deal with this by comparing the emission from an eccentric binary to the known (and more simple to calculate) emission of a circular binary emitting at the peak frequency determined from the eccentric system of interest but we allow for a power law dependence on frequency for the noise.
Note that SNR involves the ratio of amplitude to noise, whereas the peak frequency is evaluated for the emitted power. The difference is compensated by powers of n (where n labels the harmonic) and we show how to incorporate the n-dependence to get a simple formula. We show that our result is better than anticipated. We derive a fairly precise formula for (f p , e)/ (f p , e = 0) for all e ∈ [0, 1) without assuming any particular limit or using any numerical fit. We derive formulae for different power laws of the noise curve. In particular, we consider three interesting cases: . These correspond to realistic noise curves such as LISA or DECIGO at the minimum noise, lower frequency, and higher frequency. The noise curves as used at this point have essentially this form (with some small deviations for LISA). Of course real noise curves might have features that will have to be taken into account. But given the noise curves in the literature [24,25], our approximation works very well.
In this paper we clarify the original method, check in detail how well it should work in LISA and DECIGO, and extend the analysis to power law noise curves to show how accuracy can be improved even further, extending our analysis to higher frequencies as well. We show that for periodic signals, taking the sum over harmonics near peak frequency to get the signal we can to a very good approximation treat the noise curve as flat.
The plan of this paper is as follows. In Sec. 2, we present our analytical approximation based on SNR evaluated around the peak frequency, which then evolves over time. We then generalize our formulae to power-law noise curves in Sec. 3 and compare numerically the results of the approximation with those of the true noise curve for LISA and DECIGO in Sec. 4. For the latter, we allow for different redshifts (since DECIGO is sensitive to binaries at large redshifts) which simply shifts the overall test to different observed frequencies. We comment on the detection of highly eccentric binaries in Sec. 5 and conclude in Sec. 6.

LISA SNR of Eccentric BBHs
To set the stage, in this section we review the semi-quantitative arguments presented in [23], which could help to gain some intuition about the SNR for eccentric events.
We start with the time-domain formula for the SNR [22]. For a detector like LISA with noise strain S N (f ), the SNR of eccentric BBHs can be calculated as where the summation is over all harmonic components of the GWs, h n is the GW amplitude of n'th harmonic component with frequency f n = nf 0 , and f 0 is the base frequency. T O is the total observation time. For highly eccentric BBHs we need to sum over a large number of harmonics, which can be numerically challenging. However, the noise strain S N (f ) varies as a power law over much of the range of observable frequencies (except at the turnaround from one power law to another) and this observation can help to simplify the harmonic summation. In fact, it turns out to be a good approximation to pull S N out of the summation, essentially assuming constant noise. We then find that the signal would be proportional to h 2 n = h 2 ≡ h 2 c , which is simply the amplitude of GW radiation averaged over one orbit. Having such a formula allows for a much better understanding of eccentricity dependence of the signal to noise. To see how this quantity depends on eccentricity, first consider the large e limit where we keep track of = 1 − e 2 factors. Then the GW amplitude is proportional toM with M ∝ (a ) 2 , the mass quadrupole of the binary, and where a is the semi-major axis of the binary. (See the appendix for more details.) To take the time derivative, we use the fact that d/dt =ψ(d/dψ) where ψ is the true anomaly of the binary orbit on which the mass quadrupole has sinusoidal dependence -andψ ∝ (a ) −3/2 . Therefore, Notice that the four time derivatives fromM 2 introduce four factors ofψ and there is one inversė ψ from converting the time integral into an integral over the phase. So there are 3 factors ofψ ∝ (a ) −9/2 in total. In addition, the prefactor ω 0 ∝ a −3/2 . So we have the net result (a ) 4 · (a ) −9/2 · a −3/2 ∼ a −2 −1/2 . Now we rewrite in terms of f p and , so we use the formula for f p derived in [26], where m is the total mass of the binary. Below we will also use There we see if we treat the SNR as a function of f p and e, that the SNR scales with like 3/4 .
Note that this gives us only an approximate formula for the dependence on 1 − e as e → 1. As in [23], to allow for all 0 ≤ e < 1, we treat the dependence on 1 + e as an unknown to be determined numerically. As was shown in [23], a simple formula turns out to well approximate the sum over harmonics for any value of eccentricity 0 ≤ e < 1: Therefore, for binaries with little chirping during the whole observation time meaning that f p and e are relatively constant, we have, For chirping binaries, we can extend to the following generalized expression, Here f p (t) and e(t) should be calculated using Peters's equations (30) and (31), as well as (4).
In this way we can avoid the summation over GW harmonics. We check the accuracy of our approximation at LISA and DECIGO (for two different values of redshift) in Fig. 1, which shows a comparison of SNR between our formula and an explicit numerical computation summing over harmonics applicable to events at LISA and DECIGO for frequencies up to 0.1Hz and 10Hz respectively. The inclusion of the redshift affects both the original formula (2) and the approximation (7) in the same way (as expected). The only difference between the two lower panels of Fig. 1 is the shift in the initial frequency of the events as f p0 /(1 + z). We include this as a consistency check of our result. The noise curves of LISA [24] and DECIGO [25] used for this computation. UPPER RIGHT: an / num for LISA at z = 0. LOWER: an / num for DECIGO for events that occur at z = 0 (LOWER LEFT) and z = 5 (LOWER RIGHT).
Eq. (6) shows that a nonzero eccentricity will decrease the SNR for a binary with fixed peak frequency. One can again understand this suppression of SNR by noting that the eccentricity lengthens the lifetime τ for fixed f p and thus suppresses the radiation powerĖ ∼ τ −1 .
A constant noise strain turns out to be a reasonably good assumption for e < 0.9 as we can see in Fig. 1. As we will now see, the correction comes from deviations from a constant noise curve over frequency. It is nonetheless remarkable that so long as we work in terms of the measurable signal concentrated at and near the peak frequency, we can treat the noise curve as constant, even with significant power law dependence. The result is accurate to at least within a factor of two for eccentricities less than 0.9 over the entire plotted parameter range. Beyond this, for large eccentricity e ∼ 0.9, the analytical approximation and the precise formula agree less well, but still within an order of magnitude (as can be seen in the upper right panel of Fig. 5).
We see that our approximation of a constant noise curve works well in the falling part of the noise curve and the turnaround point. In the next section, we derive a formula that is even more accurate on the falling part and improves the rising part of the noise curve as well.

Precise SNR Formulae for Eccentric Binaries with Power-Law Noise
The analysis of [23] reviewed above assumed a constant noise curve over frequency. In reality, the noise curves for the two lower frequency GW detectors we are considering are projected to have a falling power law dependence (f −4 until frequencies f = 5mHz at LISA and 0.2Hz at DECIGO), where they turn around to follow an f 2 dependence. In light of this, assuming a constant noise curve might seem to be a terrible approximation. In this section we will show that the additional dependence is well under control at fixed peak frequency.
In Fig. 2 we show the (unnormalized) spectrum of eccentric binaries with f p = 0.01Hz and e = (0.1, 0.5, 0.9), along with the LISA noise curve. As e gets larger, one can clearly observe the broadening of the spectrum due to an increasing number of harmonics. We see that at larger eccentricity, the spectrum develops a secondary peak at the lower frequency end, essentially because of the 1/n factor enhancement associated with the amplitude (not the power). The key insight of the analysis below is that the power-law dependences in the noise is mathematically equivalent to taking time derivatives of the GW amplitudes. We will show that including those additional power dependencies does not heavily affect our previous result for decreasing noise S N (f ) ∼ f −4 but does bring discrepancies for rising noise S N (f ) ∼ f 2 , allowing us to correct the formula in this case as well.
We begin with the expression of the SNR in the time domain again, but this time without the harmonic decomposition, The average · · · in this expression consists of both the angular average and the time average over one orbital motion. The angular average is due to our ignorance of the direction of the orbital plane, and the time average is to make sense of the GW amplitude h c . By the standard procedure, this average can be related to an average over the mass quadrupole in the following way. (See the appendix for more details.) Here d is the distance from the binary to the detector (neglecting the redshift), ω 0 = Gm/a 3 is the orbital frequency and P = 2π/ω 0 is the orbital period. The tensor Λ ij,k is the projector onto transverse and traceless components, M ij = µx i x j is the second moment of the mass, and To compute the SNR, we now decompose the above expression into harmonics, so that the SNR becomes where h 2 c,n (t) denotes the n'th harmonic component of h 2 c (t). In general, it is the harmonic sum in this expression that makes the computation complicated. This is particularly true for eccentric binaries, because their higher harmonic amplitudes are expressed in terms of the Bessel function J n of order n, whose numerical evaluation can be rather slow.
However, the computation can be greatly simplified if S N (f n (t)) is an even integer power of the frequency, as we now show. For a space based detector such as LISA or DECIGO, we recognize three distinct regions in the noise spectrum S N (f ). First, around the minimum, the noise curve S N (f ) can be approximated by a constant up to second order effects. Second, at the low frequency part, the noise curve scales as S N (f ) ∼ f −4 . 1 Third, in the high frequency region, S N (f ) ∼ f 2 . Therefore, we shall consider these three cases in turn.
Constant noise, S N (f ) ∼ f 0 . In this case, we can take the noise strain S N (f ) out of the summation. Then we only need to compute n h 2 c,n (t) = h 2 c (t), which is nothing but the quantity in (9). Using the explicit expression for M ij andψ, we find where γ 1.1954 is the exponent that appears in the formula for the peak frequency. So we find Decreasing noise, S N (f ) ∝ f −4 . In this case we can no longer take S N (f ) out of the summation. But we observe that can be taken outside the summation. Here n p is the harmonic number of the peak harmonic. So the SNR becomes To complete the summation here, we observe that h TT n,ij ∝ cos(nω 0 t) or sin(nω 0 t). So, Therefore, So we get Increasing noise, S N (f ) ∝ f 2 . Using the reasoning similar to that above, we see that this time we need to compute Therefore, for this rising noise we have (f p , e) (f p , e = 0) S N ∝f 2 = I 2 (e).
In Fig. 3 we show the three analytical formulae for SNR, (13), (19), and (22), corresponding to three power laws S N ∼ (f 0 , f −4 , f 2 ) of the noise curve. We compare these results with explicit harmonic summation and find perfect agreement as we should since the formula is exact in this case. This figure also shows clearly why the constant-noise assumption works well even when when the noise curve scales as f −4 : This is essentially because I 0 and I −4 are numerically very similar. On the other hand, the SNR scales with e in the opposite direction when S ∼ f 2 , showing that the constant-noise assumption will not work well for the rising part of the noise curve, and an improved formula (22) was called for. We will confirm these observations in the next section by an explicit numerical check. We note that the LISA noise strain S N (f ) increases following a power law at both sides of the LISA frequency band, namely S N (f ) ∝ f 2 on the high frequency side (f > 0.1Hz) and S N (f ) ∝ f −4 at low frequencies (f < 1mHz). We can now substitute the above power law into (2) and redo the summation over harmonics. The nth harmonic introduces additional factors of powers of n for the associated time derivative and we are not left with the simple sum we had before. However, we saw that this n-dependence can be replaced by a time derivative. For a constant as it is for fixed peak frequency, any additional eccentricity dependence introduced cancels (though the overall factor can change). This is a remarkable result.

Comparing with Numerical Results
In this section we compare the SNR calculated using our approximating formulae with the numerical results obtained by directly summing over harmonics. We do this comparison for different choices of parameters and instruments with increasing generality.
First, we consider a chirping binary of black hole with m 1 = m 2 = 30M and compute its SNR collected in LISA. For this computation we assume N2A5 configuration and four years of total observation time, and assume the events arise at low redshift. 2 We calculate the ratio an / num where an is calculated using our analytical formulae, and num is from the harmonic summation. During the four-year observation, the binary will chirp (ḟ p > 0) and circularize (ė < 0 if e = 0). The rates of chirping and circularization depend on f p and e. To get a better intuition about this, we show in Fig. 4 the evolution of binaries in the (f p , e) plane over four years. To remind the reader that the actual spectra of these eccentric binaries have finite width, we show in the right panel of Fig. 4 the "width" of the spectrum of an eccentric binary. Here the "width" is defined to be the higher harmonic number at which the radiation power is half of its peak value at f p .
We note that on the decreasing part of the noise curve, the frequency changes very little due to chirping. However, at higher frequencies and, in particular, near the dip of the LISA noise curve, the chirping changes the frequency significantly. This will be important below, where the approximation works very well apart from events with high eccentricity at frequencies near but higher than the frequency at the dip of the noise curve. We consider the three analytical formulae derived in the previous section. In the two upper panels of Fig. 5 we show the result for the constant noise formula (13). We split the result into two panels, corresponding to the low-frequency/decreasing-noise and high-frequency/increasing-noise parts of the LISA, respectively. We see that the constant noise formula gives a decent estimate of the SNR for most parameter space. Except for the region with f p 0.1Hz and e 0 0.5, the constant-noise formula and the numerical result always agrees within a factor of two. As expected, the constant-noise formula works best around the minimum of the noise curve, corresponding to f p 0.01. The formula also works relatively well for the decreasing part of the noise curve. This is explained by the fact that the constant-noise formula and the decreasing-noise formula have the same limiting behavior when e → 1. The constant-noise formula does not work so well for high frequencies (f p > 0.1Hz) because the noise curve is rising here and the actual SNR has opposite asymptotic behavior as e → 1, as explicitly shown by the increasing-noise formula. Finally, there some features in the region 10 −3 Hz < f p < 10 −2 Hz and this is due to the galactic noise in the noise curve, which cannot be modeled by a simple power.
We then compare the decreasing-noise formula (19) with the numerical result in the lower-left panel of Fig. 5. We can see clear improvement compared to the constant-noise formula in the low frequency part when f p < 10 −3 Hz, where the noise curve is well approximated by f −4 so our analytical formula becomes asymptotically exact. The lack of agreement compared to the constant-noise curve for higher frequency f p 10 −2.5 Hz is due to the galactic noise, which makes the power law approximation less accurate in this region.
In the lower-right panel of Fig. 5 we show the comparison between the increasing-noise formula (22) and the numerical result. This time we see clear improvement at higher frequencies compared with the constant-noise formula. And the quality of the approximation has the opposite asymptotic behavior in the two right panels, as expected. For very large eccentricity, we see that the increasing-noise formula always tends to overestimate the SNR, because this formula assumes an increasing noise curve everywhere and thus would underestimate the noise at sufficiently low frequencies, where the actual noise curve deviates from the S N (f ) ∼ f 2 behavior. This overestimated SNR becomes more significant for higher eccentricity because a more eccentric binaries will radiate a wider spectrum extending to lower frequencies.
As can be seen from the two right panels of Fig. 5 there is a small region (f p0 ∼ 0.1Hz, e 0 > 0.8) where the analytical formula does not work as well, yielding more than a factor of two error. One needs to be careful when applying our formula to compute SNR for binaries in this region. However, here we benefit from the change of frequency due to chirping at high frequency illustrated in Fig. 4. An event moves out of the regime where the approximation works less well. Applying similar reasoning to DECIGO, the dip in noise occurs at even higher frequency where the chirping is even more significant so our approximation should apply at least as well in that frequency regime.

High Eccentricity Detection
We note a particularly interesting case when the eccentric binaries' peak frequency f p sits at the rising part of the noise curve where S N (f ) ∝ f 2 . As shown above, the SNR diverges for large eccentricity like ρ 2 ∼ (1 − e 2 ) −1 with other parameters fixed. This is due to the broadened GW spectrum at large e so that the spectrum extends to the lower frequency part where the noise is also lower. This suggests an enhanced detectability of highly eccentric binaries at the rising part of the noise curve over the naive expectation. This point was often ignored in previous studies, but can be crucial for highly eccentric binaries formed with frequencies where the detector noise curve is rising. This corresponds to f p 0.01 for LISA and f p 0.1 for DECIGO.
It is interesting to check whether we expect it to be possible to detect any eccentric binaries of this category, namely f p > 0.01 for LISA (or f p > 0.1 for DECIGO) and e > 0.9. A quick and dirty estimate for the number of observable binaries of this sort can be done by looking at LIGO events with (in principle) detectable nonzero eccentricity, namely, e > 0.01 at f p = 10Hz. It was  (19) and increasing-noise formula (22), respectively. Notice the coordinate range here extends further than in Fig. 1 From this we see that any binaries with e > 0.01 at f p = 10Hz would have e > 0.9 at the rising part of the LISA noise. Many dynamical channels predict the existence of such highly eccentric binaries. For example, [28] predicts a local merger rate 0.35Gpc −3 yr −1 for the dynamical binaries formed in dense star clusters and ∼ 9% of them would have e > 0.01 at 10Hz, corresponding to large eccentricity in the sub-Hertz region. [28] argued that around 80% percent of the KL binaries and 20% of non-KL binaries from this channel can have large eccentricity at f p 0.01Hz in LISA. There are also other channels that can produce similar highly eccentric mergers. In any case, when these mergers appear in the rising part of LISA or DECIGO noise, they can carry significantly large eccentricity, so that our method here could be quite useful.

Conclusions
A great deal of work has been devoted to detailed predictions of the wave forms for GW signals and building templates to accommodate them. However, very little attention has been paid to approximations that can be useful when predicting distributions and also when narrowing the parameter space for existing events. We have shown how to readily approximate SNR, even for highly eccentric events, in future lower frequency detectors. We summarize our result over the entire LISA frequency range in Fig. 6. Furthermore this type of analysis can sometimes highlight interesting features of the proposed detections, such as the anomalously big potential SNR for highly eccentric events.
In a future work we will follow up our analysis on relating distributions at different frequencies when the events have initially high eccentricity. In particular we can use the method outlined here to estimate the importance of intermediate frequency bands when distinguishing populations.
In this paper our main focus has been on the eccentric stellar-mass binaries. But we emphasize that our method can be directly generalized to any quasi-period sources. Notable examples include intermediate-mass black hole binaries, extreme mass ratio inspirals, and spinning compact objects. We look forward to interesting predictions for those in the future. Figure 6: The LISA SNR calculated from analytical formulae in Sec. 3 divided by the SNR calculated by harmonic summation, plotted as function of the peak frequency f 0 and the eccentricity e 0 at the initial time of observation. For this calculation we assume a binary with m 1 = m 2 = 30M at z = 0 and four years of observation time. We use (19), (13), and (22) for f p0 < 10 −2.6 Hz, 10 −2.6 Hz< f p0 < 10 −1. 3 Hz, and f p0 > 10 −1.3 Hz, respectively.

A More details
The GW amplitude is related to the mass quadrupole Q ij by Here Q ij = M ij − δ ij M kk /3 and M ij = µx i x j . The projector Λ ij,k projects a tensor to its symmetric, transverse, and traceless part, Λ ij,k (n) = P ik P j − 1 2 P ij P k , P ij (n) = δ ij − n i n j .
(Note that Q ij is traceless.) Rewriting Q in terms of M then gives Eq. (9). Then we need to know what is M and how to compute its time derivative. By rotation symmetry we can always bring the elliptical orbit into the (x, y) plane. Then, using the polar coordinates (r, ψ), the second moment of the mass M ij is M ij = µr 2 cos 2 ψ cos ψ sin ψ cos ψ sin ψ sin 2 ψ .


Knowledge management: Preserving, managing and sharing indigenous knowledge through digital library

bridge between the past and present to promote ways that are consistent with African realities.

Background
Indigenous knowledge systems (IKS) have been common in numerous African communities since time immemorial. The term indigenous knowledge (IK) is broadly defined to mean the 'creativity and dynamic ways in which people of a specific community have understood themselves concerning their natural environment' (Semali & Kincheloe 1999:3). Indigenous knowledge encompasses 'beliefs, rituals and perceptions, ways of learning, local technology, stocks of knowledge and the practices of acquiring and transmitting it' (Chambers 1991:83). In other words, it refers to a 'large body of knowledge and skills that has been built and imparted through indigenous ways that are outside the formal education system' (Lanzano 2013). Historically, IKS has been regarded as one of the indispensable assets rural people own and utilised for transforming and developing communities (Kugara 2017).
Indigenous knowledge is mainly shrouded in secrecy and the hope of it being disseminated to everyone and utilised for the public good is almost impossible. It is because of this that indigenous peoples' way of preserving, managing and sharing knowledge is argued to have contributed to their knowledge being threatened with extinction. Little or no palpable evidence exists from various studies to give knowledge that is warmly accepted to address IK management in indigenous communities. Of interest is how some writings on the phenomenon are shrouded in ignorance, chauvinism and distortion as their submissions open up for intellectual theft for indigenous peoples. The authors submit that little has been done to safeguard belittling and despising African value systems. Thus, IKS is facing extinction because of the old methods of knowledge management that are used. Mbiti (1969:194) alluded that some of the literature is full of offensive attitudes on the sensitive nature of indigenous peoples. Efforts to ignore IKS and painting it as primitive have been the foundation of some writings by some authors, thus, leaving this unique knowledge to disappear.
In Vhembe district, South Africa, the church and modern school is alleged to have played a pivotal role in the erosion of the Vhavenḓa indigenous ways of preserving, managing and sharing knowledge. The church is regarded to have destroyed the Vhavenḓa indigenous practices as they Indigenous peoples' ways of preserving, managing and sharing knowledge are argued to have contributed to their knowledge being threatened with extinction. This article addresses the need for intergenerational transmission of culturally unique knowledge for future generations through digital libraries. The speedy use of digital media in the modern, dynamic world gives valuable opportunities to facilitate the process of preserving, managing and sharing knowledge that are unique to South African communities. The use of qualitative research ensured a deep inquiry and understanding of digitalising indigenous knowledge. The modernisation theory was grounded by African indigenous ways in designing the framework. This article proposes the use of digital libraries in knowledge management so as to reach a wider audience. More so, digitalisation should be premised to ensure that the knowledge is protected against intellectual theft and allow knowledge holders and communities to benefit from such initiatives.
preached against it and called its followers 'heathens'. Songs, dances, games and other means that were once a bedrock for indigenous ways of transmitting IK were also condemned and deserted (Blacking 1967;Stayt 1931). Consequently, those who repented and denounced these so-called 'backward and uncivilised' practices were accepted into schools and secure formal jobs. By so doing, the Vhavenḓa community structures were disturbed to an extent that indigenous ways of preserving, managing and sharing knowledge were abandoned and had little or no room.
Even though customary practices play a very important role in the lives of the African people, some African governments have been reluctant to acknowledge the existence of indigenous groups within their territories. Although some of the IKS are argued to violate human rights, this cannot be a license to throw all the African value systems into the dustbin. Also, indigenous peoples' ways of preserving, managing and sharing knowledge are argued to have contributed to their knowledge being threatened with extinction because of its lagging behind of time and not conforming to the present-day ways of disseminating knowledge. Because of this, indigenous people's assets and/or ways of transmitting knowledge have become outdated and fall short of intergenerational preservation of knowledge. The available means of acquiring, preserving and sharing IK is slowly fading away thereby making it vulnerable to extinction. Resultantly, the loss of IK is highly likely to give birth to an indigenous identity crisis and a danger in recognising and understanding an 'invaluable sustainable knowledge system' in South Africa and the world at large (Human Rights and Equal Opportunity Commission 2009:221). This article can equally serve the future of African societies when it commits to preserve and learn from IK especially through digital libraries.

Literature review
The review of the literature was carried out under the following headings:

Knowledge management
The concept of IK management includes the identification, collection, codification, archiving, organising, protection, transfer, linking, application, spread and sharing of knowledge on indigenous community livelihood and environments, for sustainable development (Mabawonku 2002;Muswazi 2001). Nonaka and Toyama (2005:424) postulated that 'knowledge vision' aids to give a road map concerning the 'creation, sharing and utilisation of its knowledge'. Besides its efficacy, knowledge management at large is costly because it includes monetary, material, human and different assets for it to be successful (Davenport 1988). Kaniki and Mphahlele (2002:34) declared that most knowledge management standards can be utilised as a part of the managing of IK albeit a few issues should be tended to for IK to be viably overseen for the benefit of all. Kaniki and Mphahlele (2002:34) argued that IK was marginalised because most outsiders associated it with witchcraft, superstitions and magic. Consequently, because of this marginalisation, IK was not documented resulting in it not being viewed or recognised. As such, awareness of IK through proper management should be raised to make those who demonise it understand it before it becomes extinct and for the sake of the incoming generation (Kugara et al. 2020b).
Unlike Western knowledge that is packed and prepared for collection and housed in libraries, IK is primarily the property of communities and is shared and passed on in appropriate circumstances and as indicated by specific standards. This entails that indigenous peoples experience deep-rooted training in their specific environment, learning to subdue and work together with nature (Mumba 2002:317). Because of the nature of IK, there is a need for information experts to characterise, recognise and manage this information. Now that we live in a global village that moves with time, there is a need to get tuned in to what the current generation sensibly needs and require IK to be framed (Mumba 2002:318). In the 21st century, people have more interest in technology. Therefore, technology should be utilised to organise IK for sustainability. The current authors sought to establish whether the proposed preserving, managing and sharing strategies for IK, through a digital library, would be received by indigenous communities. Digitalising libraries for IK management entails the designing, storing and preserving of IK on digital platforms. Library as part of knowledge management is discussed in the following section.

Libraries, knowledge management and intellectual property
To start with, IK generated and acquired through libraries must be supposed to be subsequently shared: 'most libraries lack knowledge management policies and strategies to harness staff expertise for enhanced service delivery' (Swanepoel 2005:33). Furthermore, the 'organisational culture and organisational structure' are mostly not favourable for IK preserving, managing and sharing, especially with IK. Syed-Ikhsan and Rowland (2004) revealed that organisational structure in most libraries is 'protocol based' making it unfavourable and contradicts the essence of preserving, managing and sharing IK. Library use by indigenous peoples, especially digitalising, demands that users be trained to engage meaningfully in preserving, managing, sharing and accessing digitalised data.
Libraries can productively and viably design, gather, classify, sort out, safeguard and spread IK. Aside from gathering new information, libraries should also train and connect with communities. Muswazi (2001:250) and Mabawonku (2002:16) agreed that indigenous information should be preserved, managed, shared and protected by the libraries. The authors' critical opinion is that there seem to exist several potential challenges and barriers to using digital libraries to share IK. For instance, most libraries 'lack clearly defined knowledge sharing strategies' (Smith & Lumba 2008:169). In contrast, the submission that knowledge from indigenous peoples should be legally protected is tantamount to fallacy. The authors http://www.hts.org.za Open Access argue that this is problematic as there is no consensus on how to protect the intellectual property rights (IPRs) of IK practices. For this reason, fears have developed that indigenous peoples who have been in charge of development and preserving the knowledge will lose them to deceitful outsiders (Ng'etich 2005:7). This has brought about numerous arrangements of upholding the utilisation of IPRs to protect the privileges of the knowledge holders to their IK (Ng'etich 2005:8). Thus, IK is susceptible to notorious abuse and exploitation because as it stands it has not been awarded any legal protection that libraries can use.
The World Bank (2002) attested that before 1992, IK and assets were viewed as the common heritage of humankind. There were no international laws directing access to hereditary assets. Therefore, there was an expansion in the business utilisation of the knowledge and natural assets of indigenous peoples. The fast consumption of ecological assets and the need to reward the users and suppliers offered to ascend to the Convention on Biological Diversity (CBD), which out of the blue recognised the use and importance of IK. It established a system for giving access to hereditary assets and a method for fair and equitable benefit sharing (World Bank 2002:1). Indigenous knowledge has helped to fuel innovation and development in industries, ranging from agriculture and pharmaceuticals to chemicals, paper products, energy and others. For instance, Kugara et al. (2020) opined that: [F]rom time immemorial, natural disaster management in Africa has been deeply established in local communities which apply and utilise indigenous knowledge to forecast and monitor climate and other natural disasters, and build up early warnings for their benefit and future generations to come. (pp. 225-226) Intellectual property rights (IPRs) are a means for obtaining ownership for specific assets that are elusive. It, for the most part, includes the security of some type of innovation made by the human mind. This incorporates a wide assortment of creations, going from new music, novels, drugs, to computer software and products acquired from the utilisation of IK (World Bank 2002:1).
There is an emerging debate to ensure the protection of IPRs of IK practices. The United Nations Draft Declaration on the Rights of Indigenous Peoples underscores the way that indigenous peoples have the right to own and control their cultural and intellectual property relating to their science, technologies, seeds, medicines, knowledge of flora and fauna, oral traditions, designs, art and performances (United Nations 2001). Be that as it may, the conventional intellectual property law does not sufficiently cover or ensure IK and developments of indigenous peoples (Mugabe 1999). The World Intellectual Property Organization (WIPO), one of the institutions that deal with intellectual property, has just started to address these sensitive issues (Vasilev et al. 2020).

Methodology
In addition to desk-based research, data were collected qualitatively using semi-structured informal interviews and knowledge holders' experiences in two rural villages in Vhembe district, the Limpopo province of South Africa. The sampling of the study participants was both convenient and purposive. Firstly, the authors' familiarity with the study area, knowledge holders and previous study contacts simplified the data collection process. Secondly, knowledge holders were selected on the ground of their willingness to participate in the study, and all indigenous protocols of data collection and writing were meticulously followed (Cresswell & Plano Clark 2011). Key participants included traditional healthcare practitioners, elderly people (those above 60 years of age), African Studies lecturers and youth (those between the ages of 14 and 35 years). Informal focus group interview participants were sampled through a convenient sampling technique. With a convenient sampling technique, participants were involved in the study merely because they happened to be where the authors were collecting data (Etikan, Musa & Alkassim 2016). The majority of knowledge holders and participants have historically associated UNIVEN African Studies scholars and similar researchers as advocates for the development and this has effectively enhanced cooperation in data collection processes.

Findings
The study findings are presented in three themes.

The collection, preservation and dissemination of indigenous knowledge systems
The collection, preservation and dissemination of IKS were generally observed by participants to be performed in line with Sankofaism. Sankofaism is a philosophical paradigm that focuses on promoting African and indigenous ways of knowledge and skill development (Dei 2012). The African concept of Sankofa emphasises building bridges between the past and present to offer a formidable future. It provides a renaissance launch pad for the preservation of a treasured positive past, which seeks to provide the foundations that promote methodologies that are consistent with realities in Africa (Eshun 2011). It is a pivotal theme in Afrocentric education, philosophy and cultural milieus as it presents platforms for confronting mischaracterisations of African creative and pedagogical spaces (Bastos 2009 Participants were clear that the methodological process that should be followed in the collection, preservation and dissemination of IKS ought to be in line with the indigenous protocols of that community. Failure to follow that was observed to undermine the knowledge and the holders of such knowledge, thus likely to hinder its acceptance when digitised.

Indigenous ways contributed to indigenous knowledge systems being threatened with extinction
The big challenge confronting African communities is the formation of cultural intergenerational strategies to preserve, manage and share IK. This section, therefore, evaluated the efficacy of African ways of transmitting knowledge from one generation to the other. Besides this, other factors that destabilised knowledge management of African IK are colonial and postcolonial education. These two have been labelled by Whatman and Duncan (2005:120) as vital instruments that excluded positive teaching and learning of IK and maliciously labelled IK as evil and archaic. The latter is still peddled by Pentecostal churches who maintain that it is demonic. The study participants observed that features that have resulted in the loss and extinction of IK are 'institutional, structural and personal'. Most academic staff and some knowledge holders concurred that indigenous peoples' ways of preserving, managing and sharing knowledge have expressively catapulted to the erosion of this unique knowledge in the Vhembe district.
According to participants, some indigenous ways of preserving, managing and sharing IK included the use of songs, dances, games, initiation schools, storytelling, taboos and superstitions and others. One elderly observed the following regarding the usage of traditional games: 'During mahundwane indigenous knowledge was preserved, managed, and shared. Mahundwane was a game that was designed for the community children so that they be taught the norms and values of the VhaVenda people. This game was mainly and exclusively for children to learn the family setup and lifestyle management of households in the absence of parents. Elderly boys and girls who would have attended 'initiation schools' (for girls -vhusha and boys -vhutamba vhutuka) were the ones who would be in charge of the small ones to instill discipline and arrange how events unfold.' (Participant C, 2019, Vhembe District) Traditional health practitioners noticed that some IK was merely passed to certain groups of people and not everyone, underpinned to secrecy. Some of these ways were regarded as sacred as they involved the interference of the ancestors in governing its management. The youth participants were frank and outspoken that some of these practices are out of touch with time and need to be modernised as they no longer live in the past. Information that is orally transmitted is easily forgotten and subject to discrepancy each time it is passed as some aspects in the information are diminished and lost. They indicated that in each generation some things changed thereby resulting in some IK being lost. Some of these ways were argued by some parents as violating human rights, thus they were now discouraging their kids from attending. For example, Mdhluli and Kugara (2017) reports deaths of initiates that were cited as a stumbling block and that the responsible officials were dragging their feet to do away with bogus practitioners operating the 'fly by night schools'. In contrast to that, the senior citizens, although in minority, opined that the rise of misdemeanour amongst the youth created an environment that was not safe for children; hence the traditional ways of passing information for management was disrupted.
In contrast, in other societies, traditional practitioners and elderly participants observed that IK was reported to be fast phasing out, yet with serious social consequences. These participants confidently noticed that the fading away of IKS has and continues to manifest itself in numerous forms that include the rise in heinous crime rates, unemployment, mysterious missing and deaths of persons, misfortunes and diseases, moral decadence, climate change, ecosystem or environmental degradation, miscarriages and extinction of some wildlife. Even though this sparked debates with some who seemed to say they could explain the events, this group maintained that it has much to do with spirituality than what many people say. To add, it was noticed that the death of some knowledge holders and taking their knowledge to the grave without it being passed on has been a challenge with indigenous ways of preserving, managing and sharing knowledge.
Throughout history, IKS were often labelled and dismissed as archaic and 'primitive' because of the slow pace of IKS ways to conform to change. This has been mainly because of the reasons that IKS have been birthed through traditions that developed in the course of 'trial and error' and confirmed to be flexible to acclimatise with time and change (Melchias 2001:15). Even though some youth indicated some change that some come up in the music sector (incorporating indigenous values and norms), some of this 'hybrid music' was criticised by elderly participants that they are missing the point in transmitting IK as it focused on petty issues of the current youth. To clarify and shed light on the given view of Melchias, the authors beg to differ from that submission and point out that IKS were also developed by experimentations although these experiments were not documented, and the knowledge systems were legitimised and fortified under suitable institutional frameworks, culture and practices. They have been passed on to other generations (although discriminatorily) and have enabled indigenous people to survive, manage their natural resources and govern themselves.
Moreover, other participants were of the view that IKS was being examined and mercilessly judged out of context. However, these participants opined that it was a genuine and objective assessment as it sought to upgrade its ways of knowledge banking. In light of this, Lambrecht (2000) argued that the distinction between IKS and the Western knowledge system (WKS) assumes that knowledge is fixed in time, space and content whereas the creation of knowledge is a fluid process that evolves in terms of political, institutional, cultural and economic changes. The authors maintain that IKS and WKS are not fixed. Davis (2000) added that the distinction has the effect of creating hierarchies of knowledge thereby placing the so-called WKS in a dominant position and the IKS in a subordinated position. As a result, the majority of the participants agreed that to some extent IK has ceased to play its traditional role as a moral compass in the preservation, management and sharing of IK amongst indigenous livelihoods. The youth and most participants who were interviewed clamoured and advocated for modern methods that incorporate and mix African IK and scientific technology as it is more appealing to them for knowledge management.

Digital library
Indigenous knowledge has been observed to be fading away slowly in the Vhembe district and is at risk of becoming extinct. However, the authors are of the humble view that if appropriate strategies are not adopted, sooner or later this knowledge will become a history, hence the need for a digital library. The youth noticed that the advent of a digital library in preserving, managing and sharing Vhavenda IK is a way to go to restore the cultural heritage that is unknown to them. Moreover, this was labelled by the youth as an appropriate measure that has the potential of breeding new 'opportunities for development agencies, businesses, nongovernmental organisations (NGOs), and information agencies, including schools and libraries, to partner with rural communities, national governments and social entrepreneurs to create, manage and preserve knowledge and skills that are unique to communities in Limpopo. The youth were happy about this as they noticed that this has already been deemed to take place as some information centres are known for posting audio feeds and/or videos of IK to social media and communication technologies available in and around the Vhembe district. For example, Dr Rannditseni (PhD African Studies graduate) and Dr Ramavhunga (PhD African Studies graduate) have been aired out live on the Limpopo local radio engaging the members of the public with the IK they have studied. In addition, Ms Madima (PhD candidate, University of Venda -African Studies) is in the process of digitising Vhavenda indigenous games targeting the youth. It is the authors' submission that having a learned knowledge holder from the community to provide such information to the communities would illustrate the value of digital knowledge. The learned participants applauded the digital libraries as having the potency to brood an environment that permits indigenous peoples to engage and allow the osmosis of ideas about agriculture, medical care and soil conservation.
From the participants' views, an introduction and intensive expansion of social media technologies in introducing IK as a way of luring the youth to learn is a viable and foundational way of preserving, managing and sharing IK. It was unveiled in the study that there is no doubt that numerous young people are inclined to social media, thus introducing the digital library in the dissemination of IK would be a success. Moreover, these digitalised library activities should focus on the production and dissemination of standardised IK through television, telephone, Facebook, WhatsApp, radio and billboard advertisements. These digitalised IK data should appear directly to the general public, and they are expected to be fairly appropriate and particularly effective as many people today spend a great deal of their time on such media. Lor (2004) opined that 'libraries can assist with the discovery and recording of knowledge, organise it for use and promote its appreciation, including respect for the communities that produced it'. The latter submission of Lor fanned a bone of contention and forms the concept of this article and is discussed next.
The whole idea of digitalising Vhavenda IK as a way of preserving, managing and sharing for intergenerational sustenance of this unique knowledge was received with mixed feelings from other participants. Various views that seemed sensitive were captured from participants as they cautioned how and by who this should be performed. This resulted in the authors exploring the challenges and fears of knowledge holders and other participants with the idea of digitalising IK for future generations.

Challenges and fears of documenting indigenous knowledge systems
Most traditional healthcare practitioners and elderly people were passionately angered with digitalising IK for future generations as they indicated that it was going to be aired out for public consumption. They were not advocating for it as they observed that this issue had a bearing on the copyright and IPR issues that are included in the United Nations Declaration on Rights of Indigenous Peoples (United Nations Permanent Forum on Indigenous Issues 2007). Some participants suspected that the digitalising project was a western-driven project that wishes to use Vhavenda knowledge holders as tools to benefit financially from their indigenous knowledge, without sharing any benefits with the Vhavenda knowledge holders.
Participants indicated the collective nature of IK and the danger of it being shared with outsiders at the expense of an individual. In cases of intellectual property, the Economic Commission for Africa recommended that 'oral tradition and IK in African communities should be exploited in all their forms of expression, giving cognisance to the protection of intellectual property rights' (United Nations 2001:3). However, the main challenge that is a drawback in the quest to protect IK is that copyright is Eurocentric, 'placing emphasis on individuality and material', contrary to African indigenous cultural norms and knowledge (Greyling & McNulty 2011). Kuruk (2002) added that the 'historic collectivist approach to ownership of creative expressions associated with oral tradition is conceptually different from the Western copyright system, which ascribes ownership to individuals'.
Some elderly participants and the youth noticed that most library service points are located in urban areas, excluding the immense majority of the rural communities. Because the library services are geographically isolated from rural communities, it creates barriers and invisible forces that prevent rural communities from accessing library services and resources. Even wherein these library services and resources are brought to these communities, the poor infrastructure of many African libraries cannot afford the cost of digital libraries in preserving, managing and sharing IK. Besides costs, it will take a long for them to be trained to use and understand digital technology (Owiny, Mehta & Maretzki 2014).
Some participants indicated that digitalising IK would predictably result in misinformation of this unique knowledge. Mawere (2011) observed that there is a general adage that notes that 'lies, when documented, resemble truths'. Having this risk at hand, digitalising IK has been opined by elderly participants to have a high tendency of being used maliciously into cooked eye-catching inscriptions, catchy rhetorical practices and exaggerated achievements as a way of luring and winning the hearts of their targeted audiences. In that regard, Norris (1980) succinctly denounced the immorality of misinformation by opining that: [T]he advertiser should always tell the whole truth about the product he wants to sell and should judge the message not by what it says but by what the reader is most likely to think it says.
(p. 7) The digitalising of IK thus, should not manipulate, misinform or influence the prospective consumer deliberately or otherwise. Moreover, the elderly participants and traditional health care practitioners observed diversities regarding the practices should always be observed to avoid misrepresentation. This led to a discussion on how the Xhosa culture had been 'misrepresented' in a movie in South Africa and ongoing consultations with relevant stakeholders were said to be ongoing. Some Facebook comments that followed this thorny issue were as follows: Comment 1: 'If the Xhosas feel they were wrongly depicted in this film, they have aright to protest'. (sic) Comment 2: 'From an anthropological perspective, Xhosa cultural ideologies and traditions are particularly fluid and vary considerably from region to region. Also, much of it remains unrecorded (it is shared/passed on orally) so there will naturally be discrepancies. So unfortunately, it doesn't really matter how careful filmmakers were about cultural accuracy, the data will always be incorrect to at least one group or another. This scenario is a common filmmaker's trap'. (sic) Comment 3: 'Mxm…it's a movie people…do you boycott action movies because of fake bullets and explosions…This is not a documentary but a movie surely we expect some sugar and spices in it'. (sic) Because the digitalising of IK will not be for movies but for preservation, managing and sharing of IK, significant attention has to be devoted to examining the ethical problems associated with the digitalisation process. This is because some digital advertisement in South Africa has been characterised by the moral wrong of exaggerated competence.
The participants were also at pains trying to address whether digitalising of this knowledge would be possible because of the 'holistic' nature of IK. Some participants indicated that it would lose its intrinsic value as soon as it is removed from the community in which it is embedded and digitalised. Mawere (2010) alluded that: The fallacy of science is its inability to provide scientific explanations to issues of a metaphysical nature and to present alternative convincing explanations outside the canonical frameworks of scientific inquiry. (p. 13) So, whilst digitalisation strove to give convincing means of preserving, managing and sharing IK, it falls short in capturing and disseminating the unadulterated metaphysical aspect of IK. This aspect cannot be adequately answered without challenging the monopoly of science in trying to adopt it as a way of preserving, managing and sharing IK.

Conclusion
There is a unanimous growing consensus that IK is indispensable for the day-to-day running of indigenous peoples for development, governance and health. However, this unique knowledge is in danger of disappearing and becoming extinct if proper intergenerational measures are not taken to preserve, manage and share it. As such, to curb the fading away of IK a digital library is advocated for. In this article, the authors have explored the possibilities and efficacy of preserving, managing and sharing IK via digital libraries for the benefit of future generations. However, to smoothly achieve this, ethical problems associated with digitalisation process need to be meticulously observed.


Smart Topology Optimization Using Adaptive Neighborhood Simulated Annealing

Topology optimization (TO) of engineering products is an important design task to maximize performance and efficiency, which can be divided into two main categories of gradient-based and non-gradient-based methods. In recent years, significant attention has been brought to the non-gradient-based methods, mainly because they do not demand access to the derivatives of the objective functions. This property makes them well compatible to the structure of knowledge in the digital design and simulation domains, particularly in Computer Aided Design and Engineering (CAD/CAE) environments. These methods allow for the generation and evaluation of new evolutionary solutions without using the sensitivity information. In this work, a new non-gradient TO methodology using a variation of Simulated Annealing (SA) is presented. This methodology adaptively adjusts newly-generated candidates based on the history of the current solutions and uses the crystallization heuristic to smartly control the convergence of the TO problem. If the changes in the previous solutions of an element and its neighborhood improve the results, the crystallization factor increases the changes in the newly random generated solutions. Otherwise, it decreases the value of changes in the recently generated solutions. This methodology wisely improves the random exploration and convergence of the solutions in TO. In order to study the role of the various parameters in the algorithm, a variety of experiments are conducted and results are analyzed. In multiple case studies, it is shown that the final results are well comparable to the results obtained from the classic gradient-based methods. As an additional feature, a density filter is added to the algorithm to remove discontinuities and gray areas in the final solution resulting in robust outcomes in adjustable resolutions.

Introduction
Topology Optimization (TO) is a mathematical method to determine the distribution of material in a design domain in order to maximize the performance of a system [1]. The performance of each system to be optimized includes a wide range of specifications that depend on the application of each system and restrictions. The optimization process could be subjected to some constraints and the final solution must satisfy them. Depending on the TO problem physics, various methodologies have been developed in TO [2,3]. To formulate a TO problem, the design domain is discretized into Finite Elements (FE). The FE should be fine enough to describe the design domain in a reasonable resolution however it should not be too fine that it increase the computational cost without improving the results. The variables are the material density inside each element. Using the intermediate density for elements can improve some common problems: The checkerboard and the discontinuity of the design. The Solid Isotropic Material with Penalization (SIMP) method is one of the most applicable approaches to solving TO problems with less complexity. SIMP models the density of each element with a continuous value between a non-zero minimum and one. This method models the properties of material through a power law approach to penalize unfavorable intermediate densities [4].
The TO algorithm includes generating new solutions and evaluates them based on defined criteria [5]. A variety of methods have been proposed to generate new solutions and evaluate them [2]. From an optimizer point of view, TO methods can be divided into two main types: Gradient-based and non-gradient-based.
The gradient-based methods such as Optimal Criteria (OC) and the Method of Moving Asymptotes (MMA) use sensitivity information of the objective functions (that is, its first order derivatives). In order to obtain the sensitivity information, the derivation of the objective function is required. These methods have a fast convergence to the final results. The gradient-based methods are widely used for TO problems [6][7][8]. Aside from the good performance of these methods in TO, they are limited in the sensitivity information of the objective functions.
However, finding the sensitivity of the objective function is not easy [9]. In such cases, where the objective function derivative is difficult to obtain, the non-gradient-based methods are more advantageous. Non-gradient based methods work only with the value of the objective functions [10,11]. Genetic Algorithm (GA) and Simulated Annealing (SA) are two popular meta-heuristics that have been used in non-gradient-based TO. The main advantage of these non-gradient-based methods is that, unlike gradient-based methods, they do not need the sensitivity information of the objective function.
GA uses the theory of natural selection for heuristic search [12]. Despite GA being used in TO and having represented good results [13], it has high computational costs even for simple problems. The determination of some parameters is also a challenge in GA-based TO which can notably change the results. These parameters are the probability of crossover, probability of mutation, fitness scaling coefficient, population size, selection algorithm, and crossover operator [12].
SA is the other meta-heuristic for non-gradient-based TO. SA uses the concept of metal annealing for the optimization process. The SA generates new candidates randomly and accepts the one which improves the objective function. Additionally, SA also accepts the worst candidates (which do not improve the objective function) under certain conditions [14,15]. This fact allows SA to escape from local minima [16]. The Multi-Resolution Design Variable (MRDV) is used in the literature to enhance search performance [17]. Nevertheless, this approach still has high computational costs and the checkerboard pattern problem still remains. Combining SIMP with SA in structural TO is also used to obtain better results and less gray areas [18]. However, SIMP with SA is not efficient enough due to the high computational costs. Using SA in TO is almost new and needs more research to improve its convergence and computational cost which is the main objective of this work. The proposed algorithm uses the concept of SA with a crystallization factor to improve the performance of TO.
SA with crystallization heuristic controls the probability distribution of the random search [19]. A closed loop is created considering the acceptance or rejection of the new candidate. This closed loop increases the acceptance ratio of new candidates. The SA with crystallization heuristic has been successfully applied in different fields: Cutting and packing [20,21], electrical impedance tomography [22][23][24][25][26][27][28][29][30], curve fitting [31][32][33], etc. The formulation of this method and its application in TO are explained in the next section.
In this paper, the SA with crystallization heuristic has been used for TO as a nongradient-based method. The two cases of compliance minimization in the classic beam and heat transfer problem have been used to implement this method and verify the results.
Analysis of the results has also been done in order to select the input parameters such as the number of iterations and maximum temperature. The problem formulation for the optimization, SA with crystallization heuristic, and TO algorithm will be described in the next section. Then, the results from the proposed method are compared with the results from the literature.

Simulated Annealing
The SA is a probabilistic meta-heuristic inspired by the metal annealing process. SA was proposed in the early 80s based on the Metropolis algorithm [14,34,35]. The process starts with an initial solution and in each step, it generates new candidates randomly. The new candidate will be accepted if it improves the objective function or if its Boltzman probability is greater than a random number. The Boltzman probability is calculated by: where ∆E is the energy variation, which is the difference between the objective functions for the current solution and the newly generated candidate. T is the current temperature that starts from a high value and decreases in each step until reaching the frozen state. The maximum temperature should be high enough to ensure the domain exploration. Since SA can accept worst candidates, it avoids being stuck in local minima. The generation of new candidates can be performed in different ways. The crystallization heuristic proposed by Martins and Tsuzuki [36] adjusts the probability density to increase the acceptance of new solutions. The SA with crystallization heuristic determines the new candidate x new by: where x current is the current solution, ∆r i is the fixed step size associated with the i-th variable, e i is the direction of the i-th variable, and C i is the crystallization factor. Usually, The new candidate x new is determined by summing random numbers C i times (Bates distribution). The crystallization factor C i varies from 1 and a maximum value. When a new candidate is accepted, the SA considers that it is refining the solution. The SA creates a feedback action which increases exploration, resulting in the reduction of the crystallization factor (this is the negative feedback). On the other hand, when a new candidate is rejected, the SA considers that it is exploring the domain. The SA creates a feedback action which increases refining, resulting in the increase of the crystallization factor (this is the positive feedback).

Topology Optimization
Application of TO includes a wide range of problems with different physics. The problems of compliance TO in beams have been used in several research works [6,17,37]. The optimized topologies based on compliance are a good starting point for other TO problems, such as maximum stress, maximum deflection, etc. [3]. Minimizing the compliance of a structure is equal to minimizing the strain energy S which is the external work performed by the external force F with elastic deformation U from an unstressed state as shown in: For a discretized domain with N elements, the density of each element x e are the TO variables. The total strain energy can be calculated as a strain energy summation in each element with penalized design variable p. According to the SIMP method, it is calculated as: where u e is the element displacement vector and k e is the element stiffness matrix. The volume fraction constraint determines the final volume regarding the design domain. The optimization process includes the generation of new candidates and objective function evaluation. The SA with crystallization heuristic is used as an optimizer for TO. The algorithm is presented in detail in the next section.

SA in TO
In TO with SA meta-heuristic, the optimization starts by assigning an initial random solution and the initial temperature. At each step, a new candidate is generated. The new candidate is the new material density and is calculated according to (2). Then, the objective function is evaluated using FEM with the just determined material density. If the new candidate improves the objective function in comparison with the current solution, the new candidate is directly accepted. Otherwise, the Boltzman probability is calculated using (1). If the Boltzman probability is larger than a random number, the new candidate is accepted.
In both situations, if the new candidate is accepted, the crystallization factor receives negative feedback. It should be noted that the crystallization factor has 1 as the minimum value. In case the new candidate is rejected, the crystallization factor receives positive feedback. The crystallization factor has a maximum value, usually equal to 20 [38,39]. This maximum crystallization factor value ensures a minimum standard deviation value for the probability density distribution. For each temperature, the process of new candidate generation and objective function evaluation is repeated for a predefined number of iterations before going to the next temperature. The SA must reach the thermal equilibrium in the current temperature before going to the next temperature. The thermal equilibrium is reached by evaluating the objective function m times or accepting new candidates m/2 times. The value of m must be defined for each TO problem, is usually larger than 40 times the number of variables. A larger m increases the design domain exploration. On the other hand, very large m increases the computational cost, and perhaps less exploration is required. The value of m for each TO must be determined based on the number of elements (variables) and some trial runs. After reaching the thermal equilibrium, the temperature decreases by cooling factor al pha. The cooling factor is another parameter to be determined, and a value between 0.8 and 0.95 is normally used. The SA converges when the frozen state is reached, this happens when no new candidate is accepted for a specific temperature.
Since this optimization method works based on a probabilistic search, the results come with some discontinuity and gray areas. To reduce discontinuity in borders, a density filter is applied for each element. According to this filter, the density of each element is calculated by: based on the weighted density of neighboring elements, where w(x j ) is the weighting function and represents the distance of neighboring elements subtracted from the filter radius. At higher temperatures, most of the new candidates are accepted as the SA stays in an exploration phase. It is not reasonable to apply the density filter in this situation. The density filter will improve the convergence at lower temperatures, where the SA is already in the refinement phase. The transition from the exploration phase to the refinement phase is not clear. For criteria to start applying the density filter, it considered the situation where 80% or more of rejections happens to be in this transition. After reaching the convergence, a post-processing happens and the density filter is applied again to smoothen the borders and also to turn gray areas into black and white elements.

Results
In this section, three study cases are shown. The cantilever and MBB problems in 2D. These two problems allow a comparison of the obtained results with literature results. A practical 3D problem in heat conduction is also shown .

Cantilever Problem
To verify the proposed method, the classic problem of TO for the cantilever and MBB beam was selected. The cantilever beam is a design domain fixed on one end and a force applied on the other end as shown in Figure 1. The TO was applied to the cantilever beam with 60 × 30 elements with a constraint on the volume fraction. The number of elements to discrete the design domain should be high enough to give a good resolution of the final design. In this case, the number of elements was selected based on the similar cases in literature [2]. In other cases, this could be selected simply by considering the minimum required resolution of the final design which depends on the complexity of the expected final design. In order to determine the proper value for m, the optimization process for the cantilever beam problem started with m = 10. This is a very small value, the motivation was to find a value smaller than 60 × 30× 40 = 72,000, where 30 × 60 is the number of variables and 40 is the commonly used multiplier. Figure 2 shows the plot of the SA converged objective function value for the optimum solution of a 0.5-volume fraction versus the value of m.  As shown in Figure 2, the objective function (compliance) starts to converge to the optimum when m = 100. By increasing the value of m until 1000, the converged objective function value will not have any considerable changes and it is already very close to the optimum solution. Similarly, for other volume fractions, the convergence has the same behavior. Since the computational cost for this problem is not too high, m = 1000 was used in the simulations. In more complicated problems, a reasonable value of m can be determined using a similar approach.
The frozen state is defined as the temperature where no new candidate is accepted. It is necessary to determine the initial temperature value. In an attempt to stop the SA before reaching the frozen state and to define the initial temperature, the following study was executed. The number of rejected and accepted new candidates per temperature were computed. It resulted in the graph shown in Figure 3. In the SA, the number of accepted solutions must be as high as possible in the exploration phase. From Figure 3, at a temperature near 100, SA has an acceptance ratio of about 60%. The initial temperature T 0 = 100 is reasonable, and allows the SA to have enough domain exploration. The final temperature was selected as T end = 0.001 in which the number of accepted solutions is already small. The real frozen state cannot be reached, because of the checker board problem. The problem of minimizing compliance in the cantilever beam of Figure 1 is solved with the gradient-based method from the literature [6] and the method proposed herein. The results are presented in Table 1 for different volume fractions. The parameters used for the SA in the proposed method are T 0 = 100 and T end = 0.001 for initial and final temperatures, m = 1000 to reach the thermal equilibrium, and α = 0.85 for the cooling factor. The compliance of optimized cantilever beam using the proposed method agrees very well with the values from the literature. Figure 4 shows the optimized topology of cantilever beam for different volume fractions. The shapes from the method proposed herein are similar to the results from the literature.  Each variable has its crystallization factor and their value shows the crystallization of the variable as the temperature decreases. The value of the crystallization factor indicates if the variable is in an exploratory or refinement phase. Crystallization factors with higher values indicate that the variables are in the refinement phase. Crystallization factors with smaller values indicate that the variables are in the exploration phase. Four different variables were selected as shown in Figure 5a. The graph average of the crystallization factors for the four points per temperature is presented in Figure 5b. As shown in Figure 5, the variable C1 started randomly as void or solid, its crystallization factor is shown in blue and it quickly reaches the maximum value. C3 is another variable, its crystallization factor is shown in yellow. It takes more iterations to get to the frozen state. Near the temperature of 0.01, all four crystallization factors associated with the four selected variables reach the frozen state (its maximum value).

MBB Problem
The second studied problem is the Messerschmitt-Bolkow-Blohm (MBB) beam. The design domain is a simply supported beam with a vertical force in the middle. Since the design domain is symmetric, half of it considered here as showed in Figure 6. Since the design domain is completely symmetric, only half of the design domain is considered in TO. A similar procedure, as described for the previous example (see Figure 6), was performed to determine the value of m. The graph shown in Figure 7 shows the graph of the converged objective function value versus the value of m. As shown in Figure 7, with more than m = 100, the converged objective function value is near the optimum. It can be noticed from Figures 2 and 7 that the minimum required value of m depends on the number of variables, and it is almost independent of the loading and boundary conditions. The same input parameters of the cantilever beam problem are used to solve the half-MBB TO problem in Figure 6. The results for the compliance of optimized topology in different volume fractions are presented in Table 2. The results for the half-MBB problem also agrees very well with the results from the literature. The final shape for the half-MBB problem is shown in Figure 8.

Heat Conduction Problem
The results for the TO of the cantilever and half-MBB beam shows that the proposed method is as efficient as the gradient-based TO methods. It should be noted that the main significance of the proposed method is the inclusion of sensitivity information without using any gradient information. Since this method works based on a probabilistic search, the final shapes are not as smooth as the other methods. However, the use of finer mesh or the application of some post processing can improve the final results.
In order to show the usefulness of this method in practical problems, a 3D heat conduction problem was also studied. The schematic of the design domain is shown in Figure 9. In the heat transfer problem, the heat sink is at a different temperature and the heat flows throw the design domain. The formulation of the heat transfer in the FE format is similar to (3). The difference is that in the heat transfer problem U is the finite global nodal temperature vector, F is the global thermal load vector, and K is the global thermal conductivity. The thermal conductivity matrix or stiffness matrix for the heat conduction problem can be obtained by the assembly of element thermal conductivity matrices. The SIMP method gives the thermal conductivity matrix of each element by interpolating the density of each element. The design domain in Figure 9 is in 3D and all the nodes and matrices should be in 3D as well. For isotropic material, the thermal properties are equal in all directions. The objective function of this heat conduction problem to be minimized is the heat transfer compliance. According to the literature [8], a 40 × 40 element in x and y directions can define the problem efficiently. To ease the results comparison, just one element is considered in the z direction. To determine the proper number of iterations, the optimization process has been performed with different values of m and the results are shown in Figure 10. From Figure 10, the converged objective function approximated the optimum solution for values of m greater than 300. The higher the m, the slightly better the results will be, and since the computational costs are low, the value of m = 5000 was chosen. The final optimum solution values of heat transfer compliance from the literature [8] and the proposed method are presented in Table 3. The final shape of optimized topologies from the literature [8] and the proposed method are presented in Figure 11 for a 40 × 40 × 1 element design domain. As shown in Figure 11 and Table 3, the results from the proposed method and the literature for the heat transfer problem are almost the same. The value of the minimum compliance has also slightly improved.

Conclusions
A new method for TO of mechanical structures is proposed in this work. The proposed method uses the concept of SA with crystallization heuristic to determine the optimal shape. The crystallization heuristic improved the search in the generation of the new candidate. Although this algorithm is non-gradient based and derivatives of the objective functions are not required, the algorithm uses a feedback action, depending on whether the new candidate is accepted or rejected, the crystallization factor is modified accordingly. A density filter is also applied when the thermal equilibrium is reached to smoothen the borders. Applying the proposed method to the classic problem of compliance minimization in the cantilever beam and MBB beam showed good agreement with the results from the literature. Similarly, the proposed method for the problem of minimizing heat transfer compliance showed good agreement with results from the literature. It shows that this method can be used for any type of problem regardless of the existence of sensitivity information.
Since this method is new, a considerable amount of research and developments are required to improve and prepare it for commercial use. Future work could reduce the computational costs by the calculation of the objective function in different ways. This method could also be applied to multi-objective TO problems.

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

Abbreviations
The following abbreviations are used in this manuscript: