Three dimensional quantitative characterization of magnetite nanoparticles embedded in mesoporous silicon: local curvature, demagnetizing factors and magnetic Monte Carlo simulations Nanoscale

Magnetite nanoparticles embedded within the pores of a mesoporous silicon template have been characterized using electron tomography. Linear least squares optimization was used to ﬁ t an arbitrary ellipsoid to each segmented particle from the three dimensional reconstruction. It was then possible to calculate the demagnetizing factors and the direction of the shape anisotropy easy axis for every particle. The demagnetizing factors, along with the knowledge of spatial and volume distribution of the superparamagnetic nanoparticles, were used as a model for magnetic Monte Carlo simulations, yielding zero ﬁ eld cooling/ ﬁ eld cooling and magnetic hysteresis curves, which were compared to the measured ones. Additionally, the local curvature of the magnetite particles' docking site within the mesoporous silicon's surface was obtained in two di ﬀ erent ways and a comparison will be given. A new iterative semi-automatic image alignment program was written and the importance of image segmentation for a truly objective analysis is also addressed.


Introduction
Porous silicon (pSi), a versatile material achieved by nanostructuring of a silicon wafer, is employed in various elds such as in sensors, photovoltaics and energetic materials. 1 Particularly, properties like photo-2,3 and electroluminescence 4 in the visible range, its specic surface chemistry 1 and the possibility to control the pore-size without pre-structuring make it highly attractive. Its biodegradability 5 and bioactivity 6 furthermore renders this material useful for biomedical applications.
Due to the tunable size and the growth mechanism of the pores, as well as their quasi-regular arrangement, 7 pSi is used as a template for the deposition of various molecules in controlled drug delivery. 8 The embedding of iron oxide nanoparticles specically holds a lot of promise in biomedicine due to their low toxicity. Thus the combination of both materials is of fundamental interest in basic research. In the case of magnetite nanoparticle inltration within the pores, the magnetic properties can be tuned by the particle-size and their distance between each other.
Electron tomography 9,10 (ET) is a method for directly accessing the three dimensional morphology of the material with nanometer resolution. It has proven to be an invaluable tool for the characterization of nanoparticles. [11][12][13] In ET, rst a series of projection images are acquired in a transmission electron microscope (TEM) by changing the tilt angle of the sample with respect to the electron beam. Secondly the acquired tilt series needs to be accurately aligned for obtaining a high quality reconstruction. This is usually done by tracking ducial markers or by cross-correlation. In a third step, the aligned tilt series is then backprojected to form a digital three dimensional representation of the sample, usually either with a weighted backprojection (WBP) algorithm, or with a simultaneous iterative reconstruction technique (SIRT). 14,15 Finally, for any kind of quantitative analysis, the reconstructed volume is subject to a segmentation process in order to partition the volume into specic objects of interest.
Here we propose a new iterative semi-automatic alignment method, where the image under alignment is compared with its reprojection from the initial reconstructed volume, thereby showing the translation needed for a correct alignment. Also to improve the highly important segmentation step, a public library of image manipulation algorithms Insight Toolkit (ITK) 16 was employed and integrated as a dynamic library for Gatan Inc.'s Digital MicrographÔ (DM) soware. 17 To date it holds about 40 different algorithms written from the ITK library for image ltering, deconvolution, segmentation and optimization based image registration. 16 These ITK algorithms were made publicly available and have been used to process data via scripting commands in DM.
High-angle annular dark-eld scanning TEM (HAADF-STEM) imaging was used to acquire the tilt series, as it minimizes diffraction contrast, which violates the projection requirement. 18 It states that the acquired intensity/signal has to be a monotonic function of some physical property of the sample (mass/thickness in conventional TEM or Z 2 in HAADF-STEM). 19 For optimized design of a catalytic material in terms of reactivity, the knowledge of the preferred adsorption sites of the particles in terms of surface curvature is advantageous. For this reason the analysis of the local curvature (LC) distribution of the magnetite particles' adsorption site on the surface of the pSi was studied.
The demagnetizing eld generated within a particle as a response to the external magnetic eld acts to reduce the magnetic moment of the particle. This gives rise to shape anisotropy and is determined by the demagnetizing factors (DFs). For a general ellipsoid, these factors can be calculated as dened in ref. 20. Linear least squares optimization was used to t a general ellipsoid to the particles' surface mesh vertices. This allowed us to extract the area and volume of the particles, as well as the radii of the semiaxes along with their normalized directions. Therefore ET proves to be a unique tool to acquire the DF distribution of a nanoparticle assembly. The spatial and volume distribution were used as an input to magnetic Monte Carlo (MC) simulations with and without the knowledge of the DF distribution. A magnetic hysteresis curve and zero eld cooling/eld cooling (ZFC/FC) magnetization vs. T curves were simulated and compared to the measured ones from the same sample before the tomographic investigation.
The MC simulation technique with the implementation of the Metropolis 21 algorithm was chosen here, since the inuence of the nite temperature can be studied. The magnetostatic interactions of the nanoparticles in an assembly cannot be neglected due to their long range order; hence the dipolar interparticle interactions were included in the simulations. Performing simulations on magnetic nanoparticle assemblies usually requires assumptions to be made on the spatial and volume distribution of the sample. Interparticle interactions are characterized by three degrees of disorder: surface topology, volume distribution and the random distribution of the easy axes. 22 In this study, we aim to prove that ET can readily solve the rst two of these problems not forgetting the spatial distribution.
A list of abbreviations and symbols can be found in the ESI. † 2 Experimental section

Specimen preparation
The pSi templates were fabricated by anodization of a highly n-doped silicon wafer in aqueous hydrouoric acid solution. By controlling the electrochemical parameters such as current density, electrolyte concentration and bath temperature, oriented pores were grown perpendicular to the surface of the (100) silicon wafer. 23 For this work, a sample with an average pore diameter of 80 nm and a mean pore wall thickness of 40 nm was used. A thickness of about 35 mm of the porous layer has been obtained during the anodization procedure. In addition to the main pores, smaller dendritic side pores (<20 nm) were grown in the h111i direction, forming a complex three dimensional structure consisting of dendritic, oriented pores, which are separated from each other. Fig. 1(a) and (b) shows an overview of the cross-section and a high magnication scanning electron microscope (SEM) image of the investigated sample respectively before the inltration of the magnetite nanoparticles. In Fig. 1(c) the quasi-regular arrangement of the pores is shown at the surface. Fe 3 O 4 nanoparticles of $5 nm size have been inltrated into the pores, which were prepared by high temperature decomposition in the presence of an organic precursor as explained in ref. 24. To avoid agglomeration and oxidation, the particles were coated with oleic acid, and dispersed in hexane solution.
The TEM preparation of the material was accomplished using a FEI Nova 200 Nanolab Dual-Beam (FIB-SEM) Instrument, where a pillar shaped sample was milled under cryo conditions (À160 C measured at the stage) and transferred to a TEM grid as seen in Fig. 1(d) and (e). This geometry avoids selfshadowing, thickness changes during tilting and allows for higher maximum tilt angles, minimizing the missing wedge, which in turn deteriorates the resolution of reconstruction. 10 Investigation of the reconstructions of samples milled at room temperature showed that some of the magnetite nanoparticles were fully embedded inside the silicon material. Since silicon is crystalline, this implies that the ion beam was redepositing silicon in the FIB preparation process. This effect was not found in cryo prepared samples. The nal milling and polishing step was performed using low ion beam currents (1 pA). Before the TEM investigations, the sample was plasma cleaned using a Fischione Plasma Cleaner model 1020 for 2 minutes in order to minimize hydro-carbon build-up during acquisition. This procedure also effectively removed the oleic acid coating of the magnetite nanoparticles. The nal structure of the sample used in the tomographic investigations is shown in Fig. 1(f).

Measurements and equipment
From previous studies, FTIR and EDX spectra show the presence of oleic acid and magnetite nanoparticles, respectively within the pSi template. 25 Patterns of X-ray diffractograms of Fe 3 O 4 nanoparticles correspond to an inverse spinel structure with lattice parameters of 8.38Å, which is in good agreement with published magnetite patterns. 26 Magnetization measurements have been carried out with a SQUID (Cryogenics) and a VSM (Quantum Design), respectively. To achieve eld dependent magnetization curves a magnetic eld between AE6 T was applied parallel and perpendicular to the pores. Temperature dependent magnetization was measured between 4 and 300 K. 27 A 200 kV Tecnai F20 equipped with a Schottky eld-emission gun and a HAADF detector model Fischione 3000 was used to acquire a HAADF-STEM single axis tilt series from À79 to 78 with 1 tilt increments. The camera length was set to 200 mm and a spot size corresponding to 0.3 nm and a probe current of 20 pA was used. Dynamic focusing was not needed in this study because of the pillar shaped sample structure, staying in plane while tilted.

Projection-re-projection comparison alignment
The images were rst roughly aligned using the common crosscorrelation procedure followed by a method based on the center of mass as depicted in ref. 25. However, these methods were found not to be accurate enough to yield a high quality reconstruction. Hence, an iterative alignment method was implemented (realized in DM based on the ITK algorithms). Here the roughly aligned tilt series was rst reconstructed using the WBP algorithm and then re-projected to form another tilt series, which was then cross-correlated with the original to nd the correct translation. Alternatively one can calculate the maximum mutual information by exhaustive search to determine the misalignment. 28 For this procedure the term projection-re-projection comparison (PRC) seems to be appropriate.
The re-projections can be greatly enhanced by using a weighted sum of the maximum projection 29 and the usual integral projection. In maximum projection, only the brightest pixel along the projected ray is displayed, therefore making the high intensity/high contrast objects clearly visible in HAADF-STEM images, which is crucial to obtain a correct alignment. In Fig. 2 one can see zero degree re-projections taken from an initial reconstruction using Fig. 2(a) maximum projection and Fig. 2 The quality of the PRC alignment judged by a visual criterion is clearly better as seen in Fig. 2(c) and (d). Here the former image was reconstructed using SIRT with 20 iterations aer aligning the tilt series with two iterations of the usual cross-correlation procedure and one renement step based on the centre of mass. In cross-correlation, bandpass and hanning window ltering applied to a rectangular region of interest close to the tilt axis was used to optimize the process. The latter image was reconstructed approximately on the same plane, by using the same parameters. The reconstruction quality is clearly superior already aer three iterations of the PRC alignment algorithm. The global tilt axis rotation and shif have been checked between each of the three iterations.

Segmentation
The aligned tilt series was reconstructed with a SIRT algorithm over 20 iterations. For any kind of quantitative analysis concerning ET data, it is paramount to obtain an objective segmentation and to minimize the user bias. The reconstructed volume was smoothed using a modied edge-preserving curvature anisotropic diffusion lter 30 followed by a non-local means lter. To be short, the latter uses a search window to compare voxels' similarity to its neighbors, which is then used as a weighting function to be multiplied with that voxel. In the former, the strength of the smoothing parameter, i.e. the diffusion coefficient depends on the local gradient direction of the data. Diffusion takes place on voxels within a similar region rather than across boundaries, thus smoothing the data while preserving the edges. For the segmentation of the magnetite particles, the volume was rst processed using a top hat lter that takes into account the area of the objects. 31 This lter enhances the relative brightness of the nanoparticles thus allowing for a global threshold. The segmented particles were then processed by the common distance transform and watershed algorithm to label the particles and to separate the few agglomerated ones. 32 For the segmentation of the template, a script was written within the DM environment to easily dene seed points into the volume, which were used in a connected threshold algorithm. 33 This segmentation routine uses a region growing method and includes those connected voxels to the segmentation, whose intensities fall in between the user determined minimum and maximum threshold levels. Fig. 3(a) shows the generated surface from the whole segmentation. For the actual analysis only the inner part of the pSi structure was studied since the outer surface was in direct contact with the ion beam in the FIB sample preparation process. Fig. 3(b) shows a region from the central slice (xy-plane) of the reconstruction and Fig. 3(c) is the segmentation of it. The dendritic growth of the pores in the h111i direction of the silicon wafer creates a highly complex network of the pores, making it difficult to separate out the main pores quantitatively. As shown in the ESI, † the main pores investigated here were found to be separated as claimed in ref. 23. In two pores, there were no magnetite nanoparticles found at all.

Local curvature distribution
From a triangulated surface mesh (processed within the Amira 34 programming environment), the principal curvatures (k 1 , k 2 ) could be derived, opening the way for further calculations of the local curvature (LC) distribution. In particular, two approaches are proposed and discussed. In the rst method, a set of rings were dened around the particles. The rst ring goes through the branching points, which are the points shared between the particle and the template. The following rings were determined using the neighbourhood connectivity information of the triangles as illustrated in Fig. 4(d) for the rst 3 rings. The local curvature for every particle was then calculated as the mean curvature value of these points, separately for every ring. The second method involved a t of the points' three dimensional coordinates to a polynomial surface (carried out in Matlab 35 ). This method approximates the materials surface at the particle location assuming it would have formed similar to its environment without particles.
The shape index 36 can be dened as where k 1 and k 2 are the maximum and minimum principal curvatures respectively. It was chosen as a measure for the  curvature since it presents an intuitive visual interpretation of the curvature and has a singular range of values À1 # S # 1, where a positive value indicates a concave or a cap like surface and a negative value a convex or cup like surface; whereas zero depicts a near at or a saddle like surface. A similar approach to the ring method is reported in ref. 37, where the mean value of the shape index falling inside a sphere within a 1 nm radius was taken. This could lead to incorrect results if there is an adjacent surface near the particles. For both of these methods the quality of the triangulated mesh in the vicinity of sharp edges and at the transition of two different materials (the supporting template and the particle) is usually very poor. Depending on the submerging ratio r S of the surface area of the particle outside and inside the template, dened as , the mesh of the template tends to be "upward" or "downward" in the vicinity of the particle, giving positive or negative curvature values respectively, regardless of the real curvature. As a remedy, an option was set up within the Amira module to exclude rings with obviously false points; however at the expense of locality.
To obtain a more quantitative assessment on the inuence of this approach, a test image (1024 Â 1024 Â 50) was generated, where 6 pixel diameter spheres were brought into contact with a at template (Fig. 4(a) and (b)) with different submerging ratios. The local curvature of the surface should therefore be zero everywhere. A at surface was chosen here for simplicity even though the shape index is not explicitly dened. The exact number of points depends on the chosen mesh quality, on the size of the particles themselves and the submerging ratio, and differed consequently within the sample. The mesh quality was determined by a minimum edge length value, which was chosen to be non-limiting. In the test image (and in the reconstruction), the rings were separated approximately one pixel apart. The distance of the rings to the particle surface was hence directly correlated with the voxel size of the volume such that, within the pSi reconstruction, the rst ring was touching the particle (branching points), the second ring was one pixel apart, etc., amounting to 0.76 nm each. The mean value for the number of points in rings 1, 3 and 5 can be exactly given and was 38.5, 50.5 and 55 respectively. Fig. 4(c) shows the value of the shape index as a function of the submerging ratio for each ring number. The shape index for the rst ring (white squares) was always found to be positive. For the rings 2 to 4, the curvature's sign depends on the submerging ratio. Not until the 5 th ring (yellow downward triangles), the shape index was zero within a 0.05 resolution as expected. Similar results were obtained using a test template possessing positive curvature with the exemption that the curvature was always lower than the real one regardless of the ring value, as shown later.
As a comparison, we also extracted the LC distribution by tting a surface to the point coordinates with polynomial regression. 38 The mean value of the point sets' normal directions was measured to ensure the correct sign for the curvature. For a reliable polynomial t, it was crucial to include an arbitrary three dimensional rotation to the point sets, which was realized in the minimization process. Since the optimization routine did not support integer values, the optimization was performed for each polynomial degree independently in a loop, where the degree was limited to 4. Higher degrees tend to give better mean squared errors (MSE) for the distances of the original points to the tted surface, but may produce unphysical zero-crossings in the area of the particle. In order to properly include or exclude ts, threshold values were set on the MSEs. To determine the lowest order polynomial degree, individual yet minimum values for the MSE were set as thresholds, rejecting inaccurate ts or higher order polynomials. The obtained polynomial surface was then used to calculate the shape index from analytical principal curvatures and the local curvature for each particle was taken as an averaged value of a small patch around the centre point of the rst ring. In the case of the test sample in Fig. 4(c), this method yields zero curvature (black circles) for every submerging ratio. In the calculations, all the point coordinates from the second up to the h ring were included for the tting process.
In Fig. 5(a) and (b) two examples of the tting process are shown from the pSi surface. The blue line shows the mean value of the points' direction normal, that is, on which side the particle was and ensures the correct sign for the shape index aer the rotation of the point set. The coloring depicts the shape index from the t. In Fig. 5(c) one can see a triangulated mesh of one individual pore and in Fig. 5(d) a polytted surface representation of it.
The LC distributions for different ring values are shown in Fig. 6(a)-(c) acquired using the ring method. As expected, the rst ring's tted Gaussian curve consists mainly of positive values, with a mean value of 0.27 and a standard deviation of 0.09. The peak tends to settle around zero within the next rings. The distribution of the shape index for the rst method was always found to be a Gaussian curve centered at or close to zero irrespective of the number or locations of the rings used, therefore suggesting that the magnetite particles prefer near at or saddle like surfaces for adsorption, which is in accordance with the results in ref. 37.
However, the second method paints a completely different picture. There two main Gaussian curves are centered at À0. 43 and 0.51 accompanied by a smaller one at 0.07, indicating that the particles prefer convex or concave locations for adsorption. A possible explanation for the difference in the results can be found in the way the data are extracted. For the ring method, the values are closer to zero simply because of the extension of the particle, and even more, the distance of the rings from the high curvature adsorption site. This naturally explains why the two main peaks in the polytted distribution in Fig. 6(d) are both shied towards zero; forming a similar Gaussian curve as in Fig. 6(c). This is also illustrated in Fig. 4(e), where the high curvature location is correctly estimated using the polyt method. Near the particle the curvature is false because of the triangulation process and the 5 th ring gives lower curvature values than the polyt method.
In order to prove that the polyt method can accurately t high curvature locations, a test image was created (available in the ESI †). There, a particle was situated on a location with a shape index of 0.94 (Amira). With the ring method, strongly deviating results of 0.35, 0.4 and 0.5 were obtained on the rst, third and h ring respectively. In contrast to that, the polyt method yielded a value of 0.99 using the mean coordinates of the tted points. Alternatively, when taking the mean shape index value of a small patch of points (as was done in the pSi sample), i.e. points that are within a Euclidian distance of 0.5 to the mean value of the tted point coordinates, a value of 0.96 was achieved. Hence it is reasonable to conclude that the polyt method is closer to the real situation than the ring method. In Fig. 6(e) the histogram of the submerging ratios for the particles is also given.

Demagnetizing factors
A module was written to export the coordinates of the particles' point sets from Amira to Matlab environment, where an arbitrary ellipsoid was tted to the points using the linear least squares method. 39 Thresholds were used such that the surface area and volume of the particles were similar to the tted ones. The surface area of the particles was calculated in Amira as the sum of the triangle areas and the volume was calculated using the 3D convex hull of the point set. From the t it was then possible to calculate the semi-radii of the particles along with the directions of these semi-axes. The longest semi-axis can be interpreted as the shape anisotropy easy axis since the demagnetizing energy is minimized in this direction. More information on the calculation of the DFs is given in the ESI. † Essentially, the equations for the DFs can be directly applied by inserting the semi-axes (a, b and c), which are obtained from the ellipsoidal tting procedure. In Fig. 6(f) and (g) the volume distribution of the segmented particles and their ellipsoidal ts are presented respectively and the mean diameter of the particles from the former was found to be 5.9 nm. Additionally the 3D Feret ratio is shown in Fig. 6(h) as an overview of the eccentricity of the assembly. The Feret ratio was calculated by dividing the longest semi-axis with the smallest one from the ellipsoidal t.
In Fig. 7(a), the tted particles of one pore are shown, with the blue lines showing the direction of the shape anisotropy easy axis and its magnitude is the inverse of the DF of that axis. Fig. 7(b) presents the DF distribution as a polar plot, where the angular axis is the angle between the shape anisotropy easy axis and the principal z-coordinate axis. The radial axis is the value of the DF. The z-coordinate is parallel to the pores and also the direction of the optical axis in the microscope. There were no perfectly spherical particles on the red line, which draws the DF value of 1/3. However care should be taken when interpreting this plot, since for example, the closest point to the red line has semi-axes of 2.6 nm, 2.6 nm and 2.5 nm. With a 0.76 nm voxel size in the reconstruction, the 0.1 nm difference in the semiaxes becomes problematic. Also there is a clear tendency of counts in the vicinity of the zero angle with respect to the z-axis. This could be an effect of the missing wedge, which causes elongation in the same z-axis direction, although the nature of the SIRT algorithm should minimize this effect. The maximum difference of the minor diameter from the major diameter is 4.4 nm, while the average is 1.6 nm and minimum of 0.1 nm. Both the average and the maximum values are well above the error coming from the voxel size.

Monte Carlo simulations and magnetic measurements
The measured magnetic hysteresis curve and the ZFC/FC magnetization vs. T curve can be seen in Fig. 8(a) and (b) respectively. In the former, a coercive eld of 60 Oe was measured and a slight difference on the saturation magnetization between parallel and perpendicular eld to the surface can be distinguished. This is due to the difference in the spatial distribution of the particles in the parallel and perpendicular directions to the sample surface. The dipolar coupling between the particles is less in the parallel direction because the pore walls are separating the pores and the particles therein. This behavior becomes stronger as the particle radius grows. The blocking temperature T B in the latter gure was found to be 12 K, which is higher than in non-interacting assemblies of magnetite nanoparticles. The blocking temperature of 4.05 K can be calculated for the non-interacting case using K m V ¼ 25k B T B , where K m is the anisotropy energy density, V is the volume and k B is the Boltzmann constant. A mean diameter of 5.9 nm was used from Fig. 6(f) assuming spherical particles.
To obtain the magnetic conguration under an applied eld H and nite temperature T, the MC simulation, using the standard Metropolis algorithm, was employed. 40,41 More detailed information about the simulation can be found in the ESI. † At a given temperature and applied eld, the system was allowed to relax towards equilibrium for the rst 10 3 MC steps per spin, and thermal averages were calculated over the subsequent 10 4 steps. The measurements were averaged over 10 different initial conditions (random congurations of the anisotropy easy axis and initial spin orientations). The error bars were too small to be included in the gures. The mean particle volume was hVi ¼ 108.01 nm 3 according to the log-normal t in Fig. 6(f). The material values of the bulk Fe 3 O 4 , namely, the saturation magnetization M S ¼ 4.8 Â 10 5 A m À1 and the anisotropy energy density K m ¼ 1.3 Â 10 4 J m À3 were used in the simulations.
The simulations were performed in two directions for the applied eld, namely along the z-axis (perpendicular to the surface) and on the xy-plane (parallel to the surface) and by including or omitting the DFs. The xy-plane simulations were performed with the eld rotated in steps of 30 and the nal curve was taken by averaging the values over all directions. As can be seen in Fig. 8(c), where the DFs were included, there is a clear difference between the hysteresis loops taken for the two different directions of the applied eld. With the eld along the z-axis, the hysteresis loop presents a square-like shape, whereas when in plane, its shape is thinner and more inclined. Consequently, the coercive eld is much larger when the eld is perpendicular to the surface H c,perp ¼ 0.611 compared to the case of the eld parallel to the surface H c,par ¼ 0.462. In Fig. 8(d) the simulated ZFC/FC magnetization vs. T curves, performed at cooling eld H cool ¼ 50 Oe are shown. One can see that the blocking temperature is higher in the perpendicular applied eld (along the z-axis) case (T B,perp ¼ 29 K) than that of the eld applied in the xy-plane (T B,par ¼ 24.5 K). Both values were bigger than the measured one. The difference between the two directions can be explained by looking at Fig. 7(b), where a clear accumulation of counts was found to be close to the zero angle to the z-axis. This means that the demagnetizing energy in the Hamiltonian equation (see ESI †) is stronger in the z-direction.
Next, by omitting the demagnetizing (shape) energy term in the Hamiltonian, the hysteresis (Fig. 8(e)) and ZFC/FC vs. T (Fig. 8(f)) curves were calculated. The particles were assumed to be spherical and the radius was calculated from the volume distribution. The difference between the curves calculated with the two directions of the applied eld becomes smaller. The coercive eld of H c,no_demagn ¼ 0.557 and a blocking temperature of T B,no_demagn ¼ 25.5 K were found for both directions of the applied eld. Therefore the difference between the two hysteresis loops in Fig. 8(c), in this assembly model, can be attributed to the shape anisotropy. The demagnetizing energy is equivalent to a bi-directional anisotropy (shape anisotropy) with anisotropy constants K shape1 ¼ 0.5m 0 (N a À N c ) M 2 S , K shape2 ¼ 0.5m 0 (N a À N c )M 2 S and K shape1 $ K shape2 , [42][43][44] where m 0 is the permeability constant, N a is the DF of the easy axis, N c of the hard axis. Each particle in the assembly (apart from the uniaxial crystallographic anisotropy, which is randomly oriented in the assembly) has also shape anisotropy. The shape anisotropy easy axes have a narrow, non-uniform distribution around the z-axis. So adding the demagnetizing energy term in the Hamiltonian actually includes another anisotropy term (shape anisotropy), approximately directed along the z-axis. The average shape anisotropy density to the crystallographic anisotropy is K shape /K m ¼ 0.5m 0 h(N a À N c )ViM 2 S /K m hVi ¼ 1.26, meaning that the shape anisotropy plays an important role. More particles have shape anisotropy easy axis along the z-axis and this accounts for the bigger value of the coercive eld in this direction and the difference in the shape of the hysteresis loop. It explains also the increase of the blocking temperature when the cooling eld is perpendicular to the surface.
When the demagnetizing energy is omitted, the shapes of the hysteresis loops are typical of a random anisotropy assembly. The small differences between the loops along the two directions are due to the dipolar interactions between the particles in the non-symmetric shaped assembly. Considering the importance of the demagnetizing energy term in the Hamiltonian as explained above, and the better qualitative agreement of the simulated data compared to the measured ones, it is suggested that the particles are more spheroid-like. One obvious explanation for this discrepancy can be found in the missing wedge problem, which causes elongation of the particles in the z-direction making them ellipsoidal. The elongation factor in the z-direction spawning from the missing wedge is 1.15 for an averaged maximum tilt angle of 78.5 degrees. 10 If the z-coordinates of the particles are divided by the elongation factor employing a missing wedge correction (MWC) and proceed with the simulations again, the hysteresis loops as seen in Fig. 8(g) and ZFC/FC vs. T curves in Fig. 8(h) are found. Here the qualitative agreement is better than when omitting the DF's. The directions of the shape anisotropy easy axis of the particles are more randomly oriented and therefore the ZFC/FC vs. T curves taken from different directions to the surface are also more similar.
The concentration of particles is different in the vicinity of the surface and the bulk silicon and here we studied just a small fraction of the pores approximately in the middle of the sample ( Fig. 1(a)). Higher concentrations will bring exchange interactions and fanning effects. The shape of the measured FC curve at low temperatures is descending linearly, which implies weak dipolar interactions and the measured hysteresis curve does not saturate, which could be explained by a non-random assembly of ellipsoidal particles with the hard axis perpendicular to the eld.

Conclusions
Magnetic nanoparticle simulations have not been performed hitherto using ET data according to the authors' knowledge. In this work we have shown the usefulness and power of such a technique to acquire three dimensional quantitative data using fewer assumptions on the nanoparticle assemblies.
The two different methods for the distribution of the local curvature shown here are applicable to any catalytic material. The ring method is sufficient when the radius of the particle is much smaller than the radius of curvature of the surface, otherwise polynomial regression will give more reliable results provided that they are visually certied. The implementation of the ITK library commands to DM might provide an important bridge between the electron microscopy community and the biological and medical eld, the primal users of the ITK soware. Finally the introduction of the new alignment method PRC can be useful when the deposition of ducial markers like colloidal gold is not possible in practice.
The physical understanding of nature lies in the equations that are continuously tested and scrutinized yet are still not falsied. These equations are submerged in the simulations that are performed in various elds to predict and understand the response and behavior of materials and nature itself. For the eld of materials science, ET can be extremely useful as it provides nanoscale three dimensional information of real materials, which can be further used for any kind of simulations, not just limited to magnetic ones, for example the optical response of porous materials or electron transport in solar cells. The models used in these kinds of simulations for the material's morphology are oen extremely basic compared to the real acquired three dimensional data produced by ET. Therefore as a general outlook, it is the author's view that it would be highly interesting and rewarding to use ET in combination with physical simulations within various elds to produce more accurate predictions and results with fewer assumptions on the materials.