Yiran Zhao

Julie Descamps

Yoan Léger

Lionel Santinacci

Sandrine Zanna

Neso Sojic
email: sojic@u-bordeaux.fr

Gabriel Loget
email: gabriel.loget@cnrs.fr

Upconversion Photoinduced Electrochemiluminescence of Luminol-H 2 O 2 at Si/SiO x /Ni Photoanodes

Keywords: photoelectrochemistry, electrochemiluminescence, photoanodes, silicon, luminol

ABSTRACT n-Si/SiO x /Ni photoanodes are investigated for upconversion photoinduced electrochemiluminescence (PECL) employing the 3-aminophthalhydrazide (luminol)hydrogen peroxide (H 2 O 2 ) model electrochemiluminescent (ECL) system at pH 13. We show that IR excitation of the photoanode at 850 nm produces PECL emission at 440 nm resulting in an anti-Stokes shift of -410 nm at biases higher than 0 V vs Ag/AgCl. Due to the alkaline pH value of the electrolyte that prevents Ni dissolution, this upconverting system is active for a longer period (several hours) than the reported one for similar photoanodes operating at neutral pH. The PECL intensity can be controlled by the applied potential and the power density of the incident IR light. This new PECL system opens up promising perspectives for light-addressable electrodes, microscopy, and imaging.

Introduction

Electrochemiluminescence (ECL) is the generation of light in the vicinity of an electrode induced by an electrochemical reaction. This process involves a luminophore that is promoted to an excited state by an exergonic redox reaction, and which relaxes to the ground state by emitting a photon. [START_REF] Sojic | Analytical electrogenerated chemiluminescence: from fundamentals to bioassays[END_REF][START_REF] Liu | Recent advances in electrochemiluminescence[END_REF] The ECL phenomenon is widely employed for medical diagnostics and clinical assays, [START_REF] Sobhanie | Recent trends and advancements in electrochemiluminescence biosensors for human virus detection[END_REF][START_REF] Du | Single-electrode electrochemical system for the visual and high-throughput electrochemiluminescence immunoassay[END_REF][START_REF] Sakanoue | Borondoped diamond electrode outperforms the state-of-the-art electrochemiluminescence from microbeads immunoassay[END_REF][START_REF] Ma | Rational design of electrochemiluminescent devices[END_REF][START_REF] Bhaiyya | A brief review on miniaturized electrochemiluminescence devices: From fabrication to applications[END_REF][START_REF] Kang | A simple and sensitive electrochemiluminescence spectrum measurement platform and spectrum-resolved ratiometric sensor for miroRNA-141 determination[END_REF] and research is currently progressing for microscopy [START_REF] Voci | Surface-confined electrochemiluminescence microscopy of cell membranes[END_REF][START_REF] Ding | Electrochemiluminescence single-cell analysis: intensity-and imaging-based methods[END_REF][START_REF] Rebeccani | A guide inside electrochemiluminescent microscopy mechanisms for analytical performance improvement[END_REF][START_REF] Chen | Super-resolution electrogenerated chemiluminescence microscopy for single-nanocatalyst imaging[END_REF] and imaging. [START_REF] Meng | Recent advances in electrochemiluminescence imaging analysis[END_REF][START_REF] Xu | Imaging latent fingerprints by electrochemiluminescence[END_REF][START_REF] Dong | Direct imaging of singlemolecule electrochemical reactions in solution[END_REF] These applications rely mostly on the use of anodic luminophore-coreactant model systems such as tris(bipyridine)ruthenium(II) complex ([Ru(bpy) 3 ] 2+ ) with tripropylamine (TPA), [START_REF] Tokel-Takvoryan | Electrogenerated chemiluminescence. XIII. Electrochemical and electrogenerated chemiluminescence studies of ruthenium chelates[END_REF] or 3-aminophthalhydrazide (luminol) with hydrogen peroxide (H 2 O 2 ). [START_REF] Sobhanie | Recent trends and advancements in electrochemiluminescence biosensors for human virus detection[END_REF] Photoinduced ECL (PECL) is a category of ECL that is triggered at an illuminated semiconductor in depletion. In PECL, photons (λ em ) produced by the ECL reaction are emitted under irradiation with an incident light (λ exc ). [START_REF] Laser | Semiconductor electrodes. Photo-induced electrogenerated chemiluminescence and up-conversion at semiconductor electrodes[END_REF][START_REF] Luttmer | Electrogenerated Chemiluminescence: 34. Photo-Induced electrogenerated chemiluminescence and up-conversion at semiconductor electrodes[END_REF][START_REF] Zhao | Electrochemiluminescence with semiconductor (nano)materials[END_REF] In these conditions, photogenerated minority carriers (holes (h + ) in the case of an n-type semiconductor photoanode, or electrons (e -) in the case of a p-type semiconductor photocathode) [START_REF] Descamps | Anti-Stokes photoinduced electrochemiluminescence at a photocathode[END_REF] drive the ECL reaction at an onset potential which is lower than that of a non-photoactive electrode [START_REF] Zhao | Electrochemiluminescence with semiconductor (nano)materials[END_REF][START_REF] Zhao | Photoinduced electrochemiluminescence at silicon electrodes in water[END_REF][START_REF] Yu | Luminescence amplification at BiVO 4 photoanodes by photoinduced electrochemiluminescence[END_REF] and the ECL can be localized on the electrode by focusing the incident light spot. [START_REF] Laser | Semiconductor electrodes. Photo-induced electrogenerated chemiluminescence and up-conversion at semiconductor electrodes[END_REF][START_REF] Vogel | Spatiotemporal control of electrochemiluminescence guided by a visible light stimulus[END_REF][START_REF] Zhao | Metal-insulator-semiconductor anodes for ultrastable and site-selective upconversion photoinduced electrochemiluminescence[END_REF] PECL can produce Stokes (with λ em > λ exc ) [START_REF] Yu | Luminescence amplification at BiVO 4 photoanodes by photoinduced electrochemiluminescence[END_REF][START_REF] Yu | Photoinduced electrochemiluminescence at nanostructured hematite electrodes[END_REF] or anti-Stokes (i.e., upconversion type, with λ em < λ exc ) [START_REF] Laser | Semiconductor electrodes. Photo-induced electrogenerated chemiluminescence and up-conversion at semiconductor electrodes[END_REF][START_REF] Zhao | Photoinduced electrochemiluminescence at silicon electrodes in water[END_REF][START_REF] Vogel | Spatiotemporal control of electrochemiluminescence guided by a visible light stimulus[END_REF][START_REF] Zhao | Metal-insulator-semiconductor anodes for ultrastable and site-selective upconversion photoinduced electrochemiluminescence[END_REF] photon conversion depending on the experimental parameters, i.e., the semiconductor bandgap (E g ), the wavelength of the incident light (λ exc ), and the ECL system. [START_REF] Zhao | Electrochemiluminescence with semiconductor (nano)materials[END_REF] Such benefits make PECL a promising phenomenon for infrared (IR) imaging, bioanalysis, and for investigating photoelectrodes employed in solar energy conversion applications such as photoelectrochemical water splitting. [START_REF] Loget | Dispersed Ni nanoparticles stabilize silicon photoanodes for efficient and inexpensive sunlightassisted water oxidation[END_REF][START_REF] Laskowski | Nanoscale semiconductor/catalyst interfaces in photoelectrochemistry[END_REF][START_REF] Nunez | Origin of the electrical barrier in electrolessly deposited platinum nanoparticles on p-si surfaces[END_REF] Silicon (Si) is a particularly appropriate semiconductor for anti-Stokes PECL due to its abundance and its narrow bandgap (E g = 1.1 eV, indicated as a dashed line in Figure 1a) allowing to use λ exc up to the near IR. However, when employed in the oxidation regime with the aforementioned model aqueous ECL systems, Si photoanodes are subject to passivation. [START_REF] Gerischer | Photodecomposition of semiconductors thermodynamics, kinetics and application to solar cells[END_REF][START_REF] Sun | Enabling silicon for solar-fuel production[END_REF] Applying a protective coating onto Si is a strategy to solve this problem, [START_REF] Scheuermann | Design principles for maximizing photovoltage in metaloxide-protected water-splitting photoanodes[END_REF][START_REF] Hu | Amorphous TiO 2 coatings stabilize Si, GaAs, and GaP photoanodes for efficient water oxidation[END_REF] and several publications have reported outstanding properties of Si photoanodes stabilized with thin films based on Ni for the water oxidation reaction at high pH. [START_REF] Kenney | High-performance silicon photoanodes passivated with ultrathin nickel films for water oxidation[END_REF][START_REF] Sun | Stable solar-driven water oxidation to O 2 (g) Ni-oxidecoated silicon photoanodes[END_REF][START_REF] Laskowski | Junction behavior of n-Si photoanodes protected by thin Ni elucidated from dual working electrode photoelectrochemistry[END_REF][START_REF] Loget | Tailoring the photoelectrochemistry of catalytic metal-insulator-semiconductor (MIS) photoanodes by a dissolution method[END_REF][START_REF] Lee | Scalable, highly stable Si-based metal-insulatorsemiconductor photoanodes for water oxidation fabricated using thin-film reactions and electrodeposition[END_REF][START_REF] Ma | Identifying and removing the interfacial states in metal-oxide-semiconductor schottky si photoanodes for the highest fill factor[END_REF] Inspired by this research, our group reported the anti-Stokes PECL in aqueous electrolyte using metal-insulator-semiconductor (MIS) photoanodes comprising a n-type Si (n-Si) absorber substrate, a SiO x tunnel layer and a protective Ni thin film (the photoanode is denoted as n-Si/SiO x /Ni) operated with the [Ru(bpy) 3 ] 2+ -TPA model system. [START_REF] Zhao | Photoinduced electrochemiluminescence at silicon electrodes in water[END_REF] A drawback of this system is that it requires the application of relatively high overpotentials > 0.5 V under IR illumination (all potentials in this article are referred vs Ag/AgCl). [START_REF] Zhao | Photoinduced electrochemiluminescence at silicon electrodes in water[END_REF] Furthermore, to work efficiently, the [Ru(bpy) 3 ] 2+ -TPA ECL system must be operated at physiological pH, which favors the dissolution of the Ni thin film, [START_REF] Pourbaix | Atlas of electrochemical equilibria in aqueous solutions[END_REF] causing the photoanode degradation after a few tens of minutes. [START_REF] Zhao | Photoinduced electrochemiluminescence at silicon electrodes in water[END_REF] The stability issue at neutral pH has been recently addressed by replacing the Ni coating with a noble metal such as Ir or Pt. [START_REF] Zhao | Metal-insulator-semiconductor anodes for ultrastable and site-selective upconversion photoinduced electrochemiluminescence[END_REF] However, in addition to the high price of these materials, it still implies the use of relatively high (> 0.4 V) overpotentials.

Another strategy to prevent Ni dissolution is to use an alkaline ECL system, for 1a) at a low onset potential and that it retains its activity for hours.

Materials and methods

Materials and reagents

The n-type silicon wafers (1-10 Ω cm resistivity, phosphorus-doped, single side polished, 475-525 μm) (100) and the p ++ -type silicon wafers (0.001-0.005 Ω cm resistivity, borondoped, single side polished, 490-510 μm) (100) were purchased from University Wafers.

Acetone (MOS electronic grade, Erbatron from Carlo Erba) and anhydrous ethanol (RSE electronic grade, Erbatron from Carlo Erba) were used without further purification. The ultrapure water had a resistivity of 18.2 MΩ cm (Purelab Classic UV). Sulfuric acid (96%, VLSI grade Selectipur) and hydrogen peroxide (30%, VLSI, Sigma-Aldrich) were purchased from BASF and Sigma Aldrich, respectively. Luminol was purchased from Sigma-Aldrich.

Synthesis and characterization of the surfaces

The silicon wafers (p ++ -type and n-type) were degreased by sonication in acetone, ethanol, and ultrapure water for 10 min respectively. The Si surfaces were then decontaminated and oxidized in piranha solution at 105 °C for 30 min to chemically form a ~1.5 nm-thick SiO x layer, followed by rinsing with copious amounts of ultrapure water and dried under Ar flow.

The Ni thin films were deposited on the cleaned n-Si/SiO x surfaces by sputtering with a Leica EM ACE600 coating system (Ni target purity: 99.99%, Leica). The growth of the Ni film was controlled in situ using a quartz crystal microbalance, and the thickness (~2 nm) was to have a better view of most of the surface, the z-scale was restricted from 0 to 1.5 nm (the small white dots have a maximum height of ~5 nm). The chemical composition of the outermost surface was determined by X-ray photoelectron spectroscopy (XPS) using a VG ESCALAB 250 Xi spectrometer (ThermoFisher Scientific). Survey spectra were collected at a take-off angle of 90° and a pass energy of 100, using an Al K monochromated X-ray source.

Photoelectrochemical measurements

The uncoated and coated surfaces Si surfaces (1. peaks. [START_REF] Wagner | Handbook of X-ray photoelectron spectroscopy[END_REF][START_REF] Mcintyre | X-ray photoelectron studies on some oxides and hydroxides of cobalt, nickel, and copper[END_REF] The analysis of the XPS element fine spectra is described in detail in the supporting information. The analysis of the Ni 2p and the O 1s regions (Figure S1a,b) shows that the Ni film is mainly composed of NiO and Ni(OH) 2 , produced by natural oxidation under the ambient condition of this ultrathin Ni thin film. For simplification, we nevertheless refer to these electrodes as n-Si/SiO x /Ni. The black XPS spectrum was measured on a n-Si/SiO x /Ni photoanode after a long-term polarization and will be discussed later.

Results and discussion

Surface characterization

Electrochemiluminescence (ECL)

First, electrochemistry was investigated in the dark on degenerate p ++ -Si/SiO x /Ni anodes to exclude photoinduced effects. p ++ -Si, highly doped with boron atoms, was used as a substrate because it behaves as a conductor, and allows to perform electrochemical study without illumination. As schematized in the inset of Figure 1b, the ECL emission originates from the oxidation of H 2 O 2 and luminol at the Ni surface, which induces a bond-breaking mechanism where the excited 3-aminophthalate* relaxes to the ground state by emitting a photon at 440 nm. Next, PECL was investigated in the same electrolyte (10 mM luminol and 33 mM H 2 O 2 at pH 13) on n-Si/SiO x /Ni surfaces prepared with a moderately doped and photoactive n-Si wafer.

Photo-induced electrochemiluminescence (PECL)

Figure 4a shows CVs recorded in the dark on n-Si/SiO x /Ni (black curve) and under illumination at λ exc = 850 nm with a power density (P LED ) of 6.6 mW cm -2 on n-Si/SiO x (green curve) and n-Si/SiO x /Ni (blue curve). Figure 4b shows the corresponding ECL intensity profiles with the same color code. First, one can note that neither current nor ECL emission was recorded for n-Si/SiO x under illumination, showing that hole transfer is not promoted to the electrolyte phase, which is consistent with the oxidation inhibition previously observed for p ++ -Si/SiO x (Figure 3b,c green curves). n-Si/SiO x /Ni also fails in promoting charge transfer and, consequently, no ECL emission is detected in the dark (Figure 4a,b black curves).

Conversely, under IR illumination, a photocurrent is generated, which starts at -0.4 V and reaches a plateau of ~3 mA cm -2 at 0.2 V. The maximum photocurrent density is lower than that recorded with p ++ -Si/SiO x /Ni electrode, suggesting that it is limited by the photon flux at high overpotential. [START_REF] Gärtner | Depletion-layer photoeffects in semiconductors[END_REF] The simultaneously recorded PECL intensity starts around 0 V and increased considerably at ~0. Photocurrent spectroscopy was also performed on n-Si/SiO x /Ni. The blue curve of Figure 4c presents an incident photon-to-current conversion efficiency (IPCE) spectrum that was recorded on n-Si/SiO x /Ni at 0.4 V. This spectrum shows that the photoconversion efficiency is null in the near ultraviolet and starts to increase at λ exc > 360 nm. The inactivity in the UV is caused by the electrolyte absorption, plotted as a green dotted curve in the same figure. The IPCE is higher than 40% in the visible range, reaches a maximum of ~75% at 725 nm, and varies from 69 to 38% from 800 to 1010 nm. The IPCE decrease is caused by the decrease of Si absorption coefficient at high wavelengths. At 850 nm (brown dashed line), which is the excitation wavelength employed in our PECL experiments, the IPCE is 65%, corresponding to a photocurrent density value of 3 mA cm -2 for an incident power of 6.6 mW cm -2 , [START_REF] Chen | Efficiency definitions in the field of PEC[END_REF] in good agreement with the voltammetric measurements (Figure 4a, blue curve). It is also interesting to observe how P LED affects the photocurrent density (Figure 6b, black data). Again, it seems that there is a quasi-linear relationship between these two parameters, which is well in line with a light-limited photocurrent density at this potential. Accordingly, The analysis of the XPS fine element Ni 2p, O 1s, N 1s, and C 1s spectra of the electrode after polarization is described in detail in the supporting information (Figure S1a-d). This suggests that prolonged polarization in this medium induces the precipitation of an insoluble coating composed of K + (initially present in the supporting electrolyte) and organics resulting from luminol degradation products. We expect that the formation of this layer will impede photoanode activity during longer experiments.

Conclusion

In this article, we have reported the anti-Stokes PECL of the luminol-H 2 O 2 model ECL system operated in an alkaline electrolyte on n-Si/SiO x /Ni photoanodes. Here, IR excitation of the photoanode at λ em = 850 nm produces PECL emission at λ em = 440 nm resulting in an anti-Stokes shift of -410 nm, which is the largest reported so far. [START_REF] Zhao | Electrochemiluminescence with semiconductor (nano)materials[END_REF][START_REF] Zhao | Wireless anti-stokes photoinduced electrochemiluminescence at closed semiconducting bipolar electrodes[END_REF] Due to the appropriate electrolyte pH, the system can be used for a longer time than that reported for similar photoanodes at neutral pH [START_REF] Zhao | Photoinduced electrochemiluminescence at silicon electrodes in water[END_REF], as further confirmed experimentally (Figure S2), and exhibits a high IPCE. The PECL intensity can be controlled by simple physical parameters such as the applied potential and the power density of the incident IR light. According to voltammograms and OCP measurements, the photovoltage of n-Si/SiO x /Ni in this medium is ~360 mV, which is close to the 410 mV of our previous report. [START_REF] Zhao | Photoinduced electrochemiluminescence at silicon electrodes in water[END_REF] However, the emission onset is only shifted by 230 mV when comparing the photoanode with a similar but non-photoactive anode. The electrode can be used for hours and PECL decay can be conveniently restored by replacing the electrolyte. Prolonged electrolysis showed that bubble evolution caused by the degradation of the ECL reactants produces variation in the PECL intensity. Due to the high pH of the luminol-H 2 O 2 ECL electrolyte, and the fact that the photoanode system is based on Si, we think that this report is promising for future fundamental research applied to solar energy conversion. [START_REF] Loget | Dispersed Ni nanoparticles stabilize silicon photoanodes for efficient and inexpensive sunlightassisted water oxidation[END_REF][START_REF] Laskowski | Nanoscale semiconductor/catalyst interfaces in photoelectrochemistry[END_REF][START_REF] Kenney | High-performance silicon photoanodes passivated with ultrathin nickel films for water oxidation[END_REF][START_REF] Sun | Stable solar-driven water oxidation to O 2 (g) Ni-oxidecoated silicon photoanodes[END_REF][START_REF] Laskowski | Junction behavior of n-Si photoanodes protected by thin Ni elucidated from dual working electrode photoelectrochemistry[END_REF][START_REF] Loget | Tailoring the photoelectrochemistry of catalytic metal-insulator-semiconductor (MIS) photoanodes by a dissolution method[END_REF][START_REF] Lee | Scalable, highly stable Si-based metal-insulatorsemiconductor photoanodes for water oxidation fabricated using thin-film reactions and electrodeposition[END_REF][START_REF] Ma | Identifying and removing the interfacial states in metal-oxide-semiconductor schottky si photoanodes for the highest fill factor[END_REF]

Figure 1 .

 1 Figure 1. a) Normalized spectra: electrolyte absorption (10 mM luminol, 33 mM H 2 O 2 (pH = 13)) (green curve), luminol PECL emission (blue curve), infrared excitation (brown curve), and Si bandgap position (black dashed curve). The normalized PECL emission spectra of the [Ru(bpy) 3 ] 2+ -TPA system at similar photoanodes under 850 nm illumination is also shown (orange curve). Corresponding anti-Stokes shifts are shown by blue and orange arrows. b) Scheme showing the photoelectrochemical setup employed for the PECL study. Inset: charge transfer process at the electrode with the luminol-H 2 O 2 ECL system.

  measured ex-situ by atomic force microscopy (AFM). AFM pictures were acquired on an NT-MDT Ntegra microscope in semi-contact mode with FM tips (APPNANO, SPM Probe Model: FORTA-50, resonance frequency: 43-81 kHz). The images were treated and analyzed with the open-source Gwyddion software. Note that in the AFM picture shown in Figure 2a,

3 x 1 .

 1 photoelectrochemical setup is schematized in Figure1band was previously described

Figure 2 .

 2 Figure 2. a) AFM picture of p ++ -Si/SiO x /Ni. Inset: AFM height profile of the Ni thin film onto p ++ -Si/SiO x after a lift-off process. b) XPS survey spectra of p ++ -Si/SiO x (green), p ++ -Si/SiO x /Ni (blue), and n-Si/SiO x /Ni (black) after potentiostatic polarization at 0.2 V for 4.5 h under IR illumination in the ECL electrolyte.

  Figure 3b, the two ECL reactants (10 mM luminol + 33 mM H 2 O 2 ) are studied together on a Ni-free p ++ -Si/SiO x anode (green curve) and a p ++ -Si/SiO x /Ni anode (blue curve), Figure3cdisplays the corresponding ECL intensity. p ++ -Si/SiO x produces neither current nor ECL,

Figure 3 .

 3 Figure 3. a) CVs recorded in the dark on p ++ -Si/SiO x /Ni in a pH 13 aqueous solution containing 33 mM H 2 O 2 (purple curve) or 10 mM luminol (red curve). b) CVs recorded in the dark in both 33 mM H 2 O 2 and 10 mM luminol (pH = 13) on p ++ -Si/SiO x /Ni (blue curve) and p ++ -Si/SiO x (green curve). c) Corresponding ECL intensity/potential profiles. Scan rate: 50 mV s -1 .

Figure 4 .

 4 Figure 4. a) CVs recorded in the dark on n-Si/SiO x /Ni (black) and under IR illumination on n-Si/SiO x (green) and n-Si/SiO x /Ni (blue). b) Corresponding ECL intensity-potential profiles. Scan rate: 50 mV s -1 , P LED : 6.6 mW cm -2 . Inset: photograph of the surface during PECL emission. c) IPCE spectrum recorded on n-Si/SiO x /Ni (blue) at 0.4 V, overlaid with the normalized electrolyte absorption (green dotted curve, data extracted from Figure 1a). The excitation wavelength is represented by a dashed brown line. Electrolyte composition: 10 mM luminol, 33 mM H 2 O 2 (pH = 13).

  2 V. With respect to the previous situation where ECL emission was triggered by oxidation of the reactants via the transfer of e -from the liquid to the solid phase, in the present case PECL is generated via the interfacial transfer of photogenerated h + , as shown in Figure1b. It is interesting to note that the PECL curve displays a pronounced hysteresis, with the PECL emitted at the backward scan being more intense than that emitted at the forward scan. This may be partially attributed to the formation of a particularly active Ni oxide or hydroxide phase at the end of the forward scan as hysteresis was also observed for the ECL signal on p ++ -Si/SiO x /Ni, see Figure3c. The fact that the PECL hysteresis is more pronounced than the ECL could be indicative of an adaptive-type junction having physical properties evolving during operation.[START_REF] Laskowski | Nanoscale semiconductor/catalyst interfaces in photoelectrochemistry[END_REF][START_REF] Lin | Adaptive semiconductor/electrocatalyst junctions in water-splitting photoanodes[END_REF] As was the case for ECL, the PECL emission can be easily seen with the naked eye, as shown in the digital photograph presented in the inset of Figure4b. This image shows that bubbles, likely O 2 and N 2 (produced by H 2 O 2 and luminol oxidation), evolved at the solid-liquid interface during PECL.

Figure 5 .

 5 Figure 5. a) Voltammograms (thin curves) and corresponding ECL intensity/potential profiles (thick curves) recorded on p ++ -Si/SiO x /Ni in the dark (black curves) and n-Si/SiO x /Ni under IR illumination (blue curves). Scan rate: 50 mV s -1 , P LED : 6.6 mW cm -2 . b) OCP recorded as a function of time under intermittent IR illumination on n-Si/SiO x /Ni. Electrolyte composition: 10 mM luminol, 33 mM H 2 O 2 (pH 13).

Figure 5a overlays the

  Figure 5a overlays the voltammograms and the corresponding ECL intensity recorded during the forward scan at a p ++ -Si/SiO x /Ni anode (black curves) and a n-Si/SiO x /Ni photoanode

Figure 6 .

 6 Figure 6. a) Chronoamperograms (black curve) and associated averaged PECL intensity (colored disks) recorded at 0.4 V on n-Si/SiO x /Ni during the variation of P LED . b) Plots of the photocurrent density (black data) and the PECL intensity (colored data) as a function of P LED . Disks are experimental values extracted from panel a), thick black disks are j values recorded when P LED has just been applied, the thin black disks are values recorded after 1 min of irradiation at P LED = 3.3 and 5.2 mW cm -2 . The black star is the j value extracted from the CV of Figure 4a, the black curve represents the photocurrent variation expected from the experimentally determined IPCE, the dashed purple curve is a linear fit of the PECL data. c) Chronopotentiogram (black curve) and PECL intensity (purple disks) recorded on n-Si/SiO x /Ni at j = 1.8 mA cm -2 , P LED : 5.2 mW cm -2 , electrolyte flow rate: 1 mL s -1 . Electrolyte composition: 10 mM luminol, 33 mM H 2 O 2 (pH 13).

  the j values extracted from the chronoamperograms of Figure6a(Figure6b, black disks), and the CVs of Figure4a(Figure6b, black star) are relatively close to that expected from the determined IPCE value of 65%, plotted in Figure6bas a dotted black line. For the highest P LED values, j decreased after 1 min of polarization (Figure6b, thick and thin black disks for P LED = 3.3 and 5.2 mW cm -2 ), deviating from the value expected from IPCE. We attribute this decrease to the formation of bubbles at the electrode surface that reduces the active surface of the photoanode.Finally, the PECL stability of this system was tested by applying a constant photocurrent of 1.8 mA cm -2 to the n-Si/SiO x /Ni photoanode under IR illumination (P LED = 5.2 mW cm -2 ) while the potential and PECL intensity were recorded, as presented in Figure6c. The black curve shows that the potential of the photoanode varied from 50 to 400 mV over the first 1.6 h. Although its intensity varied considerably, PECL was emitted during the first hour of the experiment. After ~1 h, PECL faded and vanished. To assess if the PECL decay was caused by the degradation of the n-Si/SiO x /Ni photoanode or the electrolyte, the electrolyte was replaced at 1.6 h by a freshly prepared one, which immediately restored the initial potential value as well as the PECL intensity. After the electrolyte change, the electrode emitted PECL for ~1.5 h and a similar decay occurred, showing that PECL loss was caused by the consumption of the ECL reactants in the electrolyte. Strong variation in PECL intensity can be observed. For instance, at 0.5 h, the intensity decreases, and, at 0.8 h, increases suddenly.Shorter but similar profiles can be observed in the second part of the experiments (e.g., at 2.3 h). The fact that the PECL increases when the potential abruptly decreases suggests that such variations are caused by the release of O 2 and N 2 bubbles (as discussed previously), which alter the emission signal intensity. To get more insights into the evolution of the photoanode surface, n-Si/SiO x /Ni was analyzed by XPS after a long polarization (4.5 h at 0.2 V). As presented in Figure2b, the 2p and 2s Si peaks were no longer observable and the intensity of the Ni 2p peaks considerably decreased. Interestingly, the intensity of the C 1s peak increased and new peaks attributed to N 1s (399 eV), K 2s at (378 eV) and K 2p 3/2 (293 eV) appeared.

Acknowledgments

This work is funded by ANR (LiCORN -ANR-20-CE29-0006). Sandrine Zanna (Chimie ParisTech, Paris) and Lionel Santinacci (CNRS, CINaM, Marseille) are warmly acknowledged for their help with XPS measurements and discussions.


Guoying Dong

Yaoyao Yunlong Tang

Fiona Zhao

A 149 Line Homogenization Code for 3D Cellular Materials Written in Matlab

Keywords: Cellular material, Lattice structure, Matlab, Numerical homogenization

Cellular architectures are promising in a variety of engineering applications due to attractive material properties. Additive Manufacturing (AM) has reduced the difficulty in the fabrication of three dimensional (3D) cellular materials. In this paper, the numerical homogenization method for 3D cellular materials is provided based on a short, selfcontained Matlab code. It is an educational description that shows how the homogenized constitutive matrix is computed by a voxel model with one material to be void and another material to be solid. A voxel generation algorithm is proposed to generate the voxel model easily by the wireframe scripts of unit cell topologies. The format of the wireframe script is defined so that the topology can be customized. The homogenization code is then extended to multi-material cellular structures and thermal conductivity problems. The result of the numerical homogenization shows that different topologies exhibit anisotropic elastic properties to a different extent. It is also found that the anisotropy of cellular materials can be controlled by adjusting the combination of materials.

Introduction

Cellular materials are evolved from nature to achieve high stiffness and strength with low density. Humankind had manufactured cellular architectures such as foams and honeycombs to mimic natural cellular materials decades ago. However, they were less sophisticated than natural cellular materials due to the limitation of manufacturing techniques. Recently, with the development of Additive Manufacturing (AM), cellular materials have been fabricated with more complicated structures and novel architectures [START_REF] Schaedler | Architected Cellular Materials[END_REF]. The lattice structure, which is a type of cellular material with interconnected structs and nodes in a three-dimensional (3D) space, has been successfully manufactured by several types of AM techniques. Due to controllable properties, periodic lattice structures with different topologies and relative densities are widely used in engineering applications for high stiffness-to-weight ratio [START_REF] Queheillalt | Cellular metal lattices with hollow trusses[END_REF][START_REF] Clough | Mechanical performance of hollow tetrahedral truss cores[END_REF], energy absorption [START_REF] Schaedler | Designing metallic microlattices for energy absorber applications[END_REF], thermal management [START_REF] Wadley | Thermal applications of cellular lattice structures[END_REF], etc.

The complex structure of cellular materials not only makes it difficult for fabrications but also leads to high computational cost for the simulation of their property. Because of small features inside the material, Finite Element Analysis (FEA) requires fine mesh to model cellular materials. If the structure is constructed by lots of unit cells, the computational cost is extremely high. To reduce the complexity of simulation, homogenization methods can be used to compute the macroscopic cellular material properties. Homogenization refers to a method that can replace the composite with an equivalent material model in order to resolve the difficulty in the analysis of boundary value problem with high heterogeneities [START_REF] Hassani | A review of homogenization and topology optimization I-homogenization theory for media with periodic structure[END_REF]. In this method, the local problem is solved to obtain the homogenized material property based on unit cell. Then the overall problem is computed by assigning a solid material with the homogenized property to substitute the periodic structure.

Asymptotic homogenization has been widely used to characterize the mechanical properties of cellular materials for several decades. It assumes that each field quantity depends on two different scales, macroscopic level, and microscopic level. The basic theory has been detailed covered in some literature [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]. Bendsoe and Kikuchi [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF] used this method to get the effective elastic modulus of a unit cell. It was also implemented in a design procedure to find the optimal topology of a unit cell under a certain boundary condition. Arabnejad and Pasini [START_REF] Arabnejad | Mechanical properties of lattice materials via asymptotic homogenization and comparison with alternative homogenization methods[END_REF] investigated the mechanical properties of 6 different lattice topologies for a whole range of relative density by asymptotic homogenization. The homogenization equation is discretized and solved via FEA. This is often referred as numerical homogenization. Hassani and Hinton [START_REF] Hassani | A review of homogenization and topology optimization I-homogenization theory for media with periodic structure[END_REF][START_REF] Hassani | A review of homogenization and topology opimization II -Analytical and numerical solution of homogenization equations[END_REF] reviewed the methods of numerical homogenization to solve general boundary value problems with periodic boundary conditions. Andreassen and Andreasen [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF] summarized the formula to compute the homogenized elasticity tensor and used numerical methods to solve the two-dimensional (2D) homogenization problem by Matlab. can be computed as:

(

where is the locally varying elasticity tensor, is the volume of the unit cell, are prescribed macroscopic strain fields, while the locally varying strain fields are defined as: [START_REF] Queheillalt | Cellular metal lattices with hollow trusses[END_REF] where is the displacement fields which can be found by solving the following equation: [START_REF] Clough | Mechanical performance of hollow tetrahedral truss cores[END_REF] where is a virtual displacement field. The code provided in [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF] can solve Eq.( 3) and computes the homogenized elasticity tensor for a 2D composite materials. It can also be used to solve 2D cellular material by assigning an extremely soft second material. However, as AM has greatly relieved the manufacturing constraints in the fabrication of three-dimensional (3D) cellular materials, the interest in cellular materials is shifting from 2D to 3D. 3D cellular materials have more design freedoms so that the mechanical performance can be further optimized than 2D cellular materials. The implementation of 2D problem detailed in [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF] is not enough for the need of designers. Nevertheless, the difficulty in the analysis of cellular materials is increased by adding one dimension. The existing homogenization implementations for 3D problem mathematically complex, which hinder designers without strong mathematical backgrounds from utilizing those approaches. Therefore, it is necessary to have a simple, self-contained code that can analyze the mechanical properties of 3D cellular materials.

In this paper, an educational Matlab code is provided in Appendix A, which can compute the homogenized constitutive matrix for 3D cellular materials. It removes the barrier for designers who want to use the numerical homogenization. It also gives instructions for those who want to understand how the asymptotic homogenization method is implemented. The implementation of the homogenization code is restricted to analyzing the linear elastic property of the cellular material. The input of this code is a 3D voxel model with isotropic material properties for each element. The voxel model can be obtained by the voxelization of the geometrical model. Furthermore, to reduce the complex procedure of geometrical modeling and voxelization of 3D cellular materials, a simple and concise Matlab code provided in Appendix A can directly generate the voxel model of predefined topologies. In Section 2, the voxelization code is introduced. The format of the predefined topology will be explained, which is the input of the voxelization code. In Section 3, the implementation of numerical homogenization will be illustrated. Section 4 discusses the result obtained from the homogenization code. In Section 5, some extension of the Matlab code will be introduced to broaden the application. Finally, this paper is wrapped up in Section 6 with some future work.

2.Voxel Generation

Numerical homogenization requires a finite element model of the lattice structure. That means a geometrical model should be generated at first. Then, the geometrical model is discretized by 3D solid elements. These procedures are tedious and time-consuming. Another way is to use a voxel model to represent the lattice structure. This still requires the generation of the geometrical model. To facilitate the procedure of getting the homogenization model of lattice structures, a voxelization algorithm based on the lattice wireframe has been proposed. In this algorithm, the voxels are directly generated from the lattice wireframe, which skips the generation of the geometrical model. This voxel model can be used as an input to the numerical homogenization code.

Read the information of lattice wireframe (lines 50-72)

To generate the voxel model, the information of the lattice topology wireframe needs to be pre-defined. The wireframe contains two parts, GRID, and STRUT. GRID defines the coordinates of all the node in the wireframe. Its format is shown in Table A1 in Appendix C, where 'ID' is the grid index; 'x' is the x-coordinate; 'y' is the y-coordinate and 'z' is the z-coordinate. STRUT defines the start grid index and end grid index. Its format is shown in Table A2 in Appendix A, where 'ID' is the strut index; 'Start' is the grid index of the start node; 'End' is the grid index of the end node. These grid indices refer to the 'ID' in the GRID. Each field of GRID and STRUT take 8 spaces. The wireframe script can be downloaded from the link in Appendix A. Furthermore, users can define new wireframe script based on the format shown in Table A1 and Table A2 in Appendix A. The visualization of these topologies is shown in Figure 1.

Generate the voxel in the design space (lines 7-21)

The design space in this research is dimensionless, which is inside a cube with a unit length. By dividing the unit length with the number of voxel on each axis, the size of each voxel can be obtained. To store the information of each voxel, the index and the center coordinates are calculated and saved. The index is used to identify the position of the voxel in the 3D logical matrix. As shown in Figure 2(a), the voxels are all generated in the first octant of the 3D Cartesian coordinate and they start from the origin. Figure 2(b) shows voxel in each layer along the z-axis. The gray voxel means it has material while the white color means there is no material in this voxel. The 3D logical matrix that represents this voxel model is shown in Figure 2(c). The coordinates measure the center position of each voxel which is used to calculate the distance between the voxel and the lattice strut. Once the number of voxels along each axis is determined, the size, index, and position of each voxel can be generated.

Determine if the voxel has material (lines 23-46)

After the generation of voxels, the next step is to determine if the voxel has material. The criterion is to measure the minimum distance from the center of the voxel to each line segment of the wireframe. If the distance is less than the radius of the lattice strut, that means this voxel is inside the lattice. Therefore, the value of this voxel will be set to '1', which means it has material. Otherwise, the value of this voxel will be set to '0', which means it is empty. When calculating the minimum distance, two cases should be considered as shown in Figure 3. Two planes, Plane A and Plane B are defined to distinguish Case (1) and Case [START_REF] Queheillalt | Cellular metal lattices with hollow trusses[END_REF]. Plane A and Plane B are located at the start point and the end point of the line segment, respectively. Both planes are perpendicular to the line segment. In Case (1), the center of the voxel is located inside the area between Plane A and Plane B. In this case, the minimum distance is from the center point to the line. However, if the center of the voxel is located outside the area between Plane A and Plane B, which is Case (2), the minimum distance is from the center point to the start point or the end point of the line segment. If the center point is exactly on Plane A or Plane B, it is included in Case (2). To decide whether a voxel is in Case [START_REF] Schaedler | Architected Cellular Materials[END_REF] or Case (2), two angles, and are defined in Figure 3. If both and are acute angles, this voxel is in Case [START_REF] Schaedler | Architected Cellular Materials[END_REF]. Otherwise, it is in Case [START_REF] Queheillalt | Cellular metal lattices with hollow trusses[END_REF]. Therefore, the first step to determine whether a voxel has material or not is to calculate and . Then, the minimum distance between the voxel center and the ling segment is determined. The last step is to compare the minimum distance with the radius of lattice strut. If the minimum distance is smaller than the radius, the index of this voxel will be set to '1'. Otherwise, it is set to '0'. After calculating all the voxels, the 3D logical matrix represents the voxel model of lattice structures is obtained. Examples of voxel models for seven topologies of lattice structures are shown in Figure 4. These models can be directly used for the proposed homogenization method in Matlab.

Line segment of lattice wireframe

3.Matlab Implementation for 3D Homogenization

In the implementation of numerical homogenization by Matlab, there are three inputs (lx), (ly) and (lz) indicate the length of the unit cell along x, y, and z-axis, respectively. The input arguments (lambda) and (mu) are the Lame's first and second parameter for the solid material which the lattice structure is made of. The last input argument (voxel) is the material indicator matrix, which specifies whether a finite element contains materials ( ) or not ( ). This matrix can be obtained by the method proposed in Section 2. It requires the finite elements with materials to be connected. If they are not connected, a two materials model with a very soft second material should be used, which will be introduced in Section 5.

Stiffness matrix and load vectors of hexahedron (lines 99-149)

The type of element used in this homogenization is an 8-node hexahedron element [START_REF] Chandrupatla | Introduction to finite elements in engineering[END_REF].

The local node number of hexahedron element is shown in Figure 5. The natural coordinates for it are called , , and . These coordinates are also shown in Figure 5. The shape function can be summarized as:

where , and are the natural coordinates of the node. The element stiffness matrix is given by the standard formula: [START_REF] Wadley | Thermal applications of cellular lattice structures[END_REF] where B is the strain-displacement matrix, C is the constitutive matrix for the element, which is constant over the element. The material is assumed to be isotropic over each element. The constitutive matrix can be expressed by Lamé's parameters as: [START_REF] Hassani | A review of homogenization and topology optimization I-homogenization theory for media with periodic structure[END_REF] where and are the first parameter and second parameter of Lamé's constants. They can be computed by the following formula [START_REF] Lautrup | Physics of continuous matter[END_REF]: [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] where E is Young's modulus and is the Poisson's ratio. By using Lamé's parameters, the element stiffness matrix can be split into two parts as [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF] which will simplify the extension to two or more materials. However, it should be noted that Eq. ( 8) is restricted to isotropic materials. The load vector is calculated by the macroscopic volumetric straining [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF] where the macroscopic strains are [START_REF] Arabnejad | Mechanical properties of lattice materials via asymptotic homogenization and comparison with alternative homogenization methods[END_REF] Similarly, the load vector can also be split into two parts as [START_REF] Hassani | A review of homogenization and topology opimization II -Analytical and numerical solution of homogenization equations[END_REF] The split parts of element stiffness matrix and load vector are calculated by calling the function 'hexahedron', which will output , , and .

Periodic boundary conditions (lines 22-39)

The structure of the mesh to model the unit cell is illustrated in Figure 6. There are 8 elements in this example. The ID of each element is shown in Figure 6(a). The nodes can be divided into three layers. Firstly, the node number without the periodic boundary condition for each node are shown in Figure 6(b). As for each element, the 24 degrees of freedom are saved in 'edof'. Each row of 'edof' is for one element and the bold number represent the element Id. The sequences of degrees of freedoms in each row of 'edof' are arranged according to the local degrees of freedoms. To demonstrate it, the comparison between global node numbers and local node numbers for the first element is shown in Figure 6(c). Secondly, lines 24-32 generate the periodic boundary condition by mirroring the back, left and top border in sequence and the periodic degrees of freedom are shown in Figure 6(d). The relationship between the non-periodic and periodic degrees of freedom is obtained in lines 34-37 and stored in 'dofVector'. Finally, line 38 updates the degrees of freedom of each element from non-periodic to the periodic condition.

Global stiffness matrix and load vector (lines 42-55)

Before assembly of the global stiffness matrix, the indices of the non-zero entries in the global stiffness matrix are created in two vectors by line 42 and 43. These two vectors are used to denote the position of the entry in the global stiffness matrix. For single material lattice structures, there are two types of elements. One type is the solid part that will be assigned to material properties. The other type is void in the lattice structure, where the Lamé's parameters are zero. Line 45 creates the matrices with material properties and for each element by the indicator matrix 'voxel'. These matrices are then used to compute the stiffness matrix for each element by multiplying the material properties and with and in line 47. After obtaining the stiffness matrix for each element, the sparse global stiffness matrix is assembled by indices vectors and element stiffness matrix in line 48. Due to rounding errors, the global stiffness matrix may not be exactly symmetric. Line 49 can eliminate the error to make the matrix symmetric so that the speed of the solver can be accelerated. For 3D problems, the computational cost could be a problem. Therefore, the symmetry of global stiffness matrix is important in solving the large system of linear equations.

In the assembly of load vectors, the indices vectors 'iF' and 'jF' are generated in lines 51-52 to denote the row and column of non-zero entries of the global load vector. For the element with solid materials, the local load vector is calculated by multiplying the material properties and with and in line 54. For the elements that are void, the load vectors are zero. Finally, the global load vector is assembled in line 55.

Solution by PCG method (lines 59-66)

Compared with two-dimensional problems, the 3D problem considerably increases the computational cost. Using in Matlab to solve global equation cost a lot of memory if the number of elements is large. Preconditioned conjugate gradients (PCG) method is used to solve the system of linear equations. Even though it not as fast as using in some cases, it can save lots of memory. The comparison between these two methods are shown in Table 1. It is found that both methods have advantages in different cases. In this code, PCG is the default method. However, is also an optional method which is deactivated in lines 65-66. If users need to use the direct method, they can simply activate lines 65-66 and deactivated lines 60-64. It should be noted that when using PCG method, the matrix should be symmetric and positive definite. The global stiffness matrix satisfies these requirements so that it can be efficiently solved by PCG method. To accelerate the speed of computation, only the degrees of freedom of element with solid materials are activated in line 58. The void element without material properties will not be considered in the computation. But this procedure requires these active elements to be connected. To accelerate the convergence of PCG method, a proper preconditioner should be obtained at first. The function 'ichol' with only one input argument is used to construct an incomplete Cholesky factorization with zero fill. The output can be used as a preconditioner in PCG method. Furthermore, the PCG method should specify the tolerance and the maximum number of iteration. If the relative residual is less than the tolerance or if the number of iteration exceeds the maximum value, the PCG method will stop and return a result with the minimum relative residual. In this paper, the default tolerance is and the maximum number of iteration is 300. These two parameters can be adjusted by users. For instance, a lower tolerance and a larger number of iteration will return a more accurate result but it takes more time to run the PCG method. In lines 61-64, PCG is run six times to get the displacement field under six loading conditions.

Homogenization (lines 69-96)

The homogenization procedure is similar to the procedure for two-dimensional problems proposed by Andreassen and Andreasen [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF]. Firstly, the element displacements corresponding to unit strain cases are found in lines 69-82. Elements' nodal displacement are solved under the uniform strains in Eq. ( 3). The difference between the 2D and 3D problem is the number of constrained degrees of freedom. The 2D problem only needs to constrain 3 degrees of freedom, while the 3D problem needs to constrain 6 degrees of freedom so that the element stiffness matrix is not singular. Line 72 explains which degrees of freedom are constrained to calculate the elementary displacements. Because all the elements are equivalent, the resulting displacements are the same for each element.

The next step is to use the element displacements and the global displacement field to calculate the homogenized constitutive matrix . The entries in can be calculated by the following equation: [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF] where is the element displacements corresponding to the th unit strain in Eq. ( 10).

is the displacement field obtained from the global stiffness equation. is the total volume of the unit cell of the lattice structure. The summation is computed in line 85-96. The part and part are computed separately and these two parts are multiplied with and . Then, the result for all the element are summed together to get . After the iteration for all six unit strains, the homogenized constitutive matrix is obtained.

Results and Discussion

Homogenized constitutive matrix

In this sub-section, an example is used to illustrate how to use the proposed method to get the homogenized property of lattice structures. To run the homogenization code, the voxel model has to be generated first. The lattice structure with the cubic topology is used as an example. The strut radius of the lattice structure is set to be 0.1; the path of the strut file is saved in variable 'address'; the resolution is set to be 50. Run the file 'GenerateVoxel.m' that can be downloaded in Appendix A as:

[voxel,Density] = GenerateVoxel(50,address,0.1);

where 'voxel' is a logical matrix that indicates the element in the lattice structure has materials or not. 'Density' estimates the relative density of this lattice structure, which is 8.44% in this case. The next step is to determine material property. It is assumed that the Young's modulus of the material is 200GPa and the Poisson's is 0.3. From Eq.( 7), the Lame's first and second parameters are calculated as 115.4 and 76.9, respectively. The size of the unit cell is 1 in x, y, and z-direction. Run the file 'homo3D.m' that can be downloaded in Appendix A as: The homogenized constitutive matrix is 0.14 0 0 0 0 0 0 0.14 0 0 0 0 0 0 0.14 0 0 0 0 0 0 7.06 0.37 0.37 0 0 0 0.37 7.06 0.37 0 0 0 0.37 0.37 7.06 [START_REF] Chandrupatla | Introduction to finite elements in engineering[END_REF] From the homogenized constitutive matrix, it can be concluded that this lattice structure has orthotropic material properties. Because of symmetry, the property along x, y and z directions are the same. To get the Young's modulus along the axial direction, the inverse of the constitutive matrix should be calculated, which is the homogenized compliance matrix . In the compliance matrix, it has: [START_REF] Lautrup | Physics of continuous matter[END_REF] where , and are the Young's modulus along x, y, and z-direction. From Eq. ( 14), it can be computed that . Because the material property is anisotropic, the Young's modulus is different along different orientations. To get the Young's modulus in other directions, the constitutive matrix has to be transformed. Firstly, the constitutive matrix is rewritten as a fourth-order stiffness tensor. Then the tensor is transformed with a rotation matrix. The transformed tensor can be used to get elastic modulus along the new direction. Through this procedure, directional properties of a cellular can be obtained.

CH = homo3D(1,
The proposed method can also be used to homogenize the periodic cellular materials with the randomized unit cell as shown in Figure 7. The length of the randomized unit cell is 2 in x-axis and 1 in y, z-axis. Therefore, the input of the unit cell length in the code should correspond to designed length. The dimension of the voxel model of the unit cell is as shown in Figure 7(c). Run the file 'homo3D.m' that can be downloaded in Appendix A as:

CH = homo3D(2,1,1,115.4,79.6,voxel);
The result of this type of cellular material is .

Effective Young's modulus surface

To further illustrate how the homogenized constitutive matrix can be used to predict the elastic properties of lattice structure in different directions. Seven types of topologies listed in Figure 1 have been analyzed by the proposed Matlab code. The bulk material properties are the same as those used in Section 4.1. For each topology, the voxel model of lattice structure is generated with three different relative density: 10%, 30% and 50%. Because the total strut lengths in different topologies are different, the strut radius need to be determined case by case. Table 2 shows the strut radius for each topology with different relative densities. The strut radius is dimensionless and the length of the unit cell is regarded as 1. It can be found that the Cubic and X-shape topologies exhibit the higher anisotropy than other topologies. Because for isotropic materials, the Young's modulus surface should be a spherical surface that means the Young's modulus is identical in all directions. However, the Cubic topology has much higher Young's modulus in axial directions while the X-shape topology has higher Young's in diagonal directions. When the relative density is 10%, the strength of these two topologies in weak direction is close to zero, which means they are very soft in those directions (diagonal direction for Cubic and axial direction for X-shape). However, with the increase of relative density, the anisotropy in these two topologies is decreasing. For instance, the maximum Young's modulus is more than 10 times larger than the minimum Young's modulus of the Cubic topology when relative density is 10%. But the ratio is reduced to around 2:1 when the relative density increases to 50%. The X-shape topology follows the similar trend. Even though they exhibit high anisotropy, the maximum Young's modulus in these two topologies is higher than other topologies. Therefore, if the principal stress is in a certain direction, these two topologies might be preferred. But if the loading case is varying and the stress field is complex, these two topologies with a low relative density should be avoided.

The material property of Tesseract topology also exhibits high anisotropy when the relative density is low. But when the relative density reaches to 50%, it is close to isotropic materials. As for Cubic-center, Face-center, Octet and Vintiles topologies, they are less anisotropic compared to other topologies. For example, when the relative density is 10%, the Face-center topology has a maximum Young's modulus around 5.5 GPa and a minimum Young's modulus around 4.8 GPa. It is suitable for complex loading cases because the stiffness in all the directions is higher than other topologies. Furthermore, the anisotropy of these four types of topologies does not reduce much with the increase of relative density. It is also found that these topologies have a lower maximum Young's modulus than Cubic and X-shape topologies. However, the Tesseract topology has a lowest maximum Young's modulus in all topologies. Therefore, if the design objective is a high stiffness, it should be avoided to use Tesseract topology. Another interesting finding in Figure 10 is that these topologies can be classified into two groups according to the direction of their maximum Young's modulus. The Cubic, Face-center, Tesseract, and Vintiles topologies have a higher stiffness along axial directions. On the contrary, the X-shape, Cubic-center and Octet topologies have a higher stiffness along diagonal directions. As a result, designers can choose the topology based on the direction of the principal stress if a higher stiffness is required. The result from the homogenization code can be used as a guide for designers to make decisions.

Convergence of the result

It is known that a higher resolution will result in more computational costs. Appropriate voxel size must be chosen to ensure accurate results with allowable computational costs.

In this subsection, the relationship between the resolution of the voxel model and the convergence of the result will be discussed. The Cubic-center topology is used in this investigation. The resolution is set from 20 to 100 with 5 as an increment. The resolution means the number of voxels in the voxel model along x, y, and z-direction. The axial Young's modulus is selected to represent the result. The bulk material property is also the same as in Section 4.1. The homogenized constitutive matrix is computed for two relative densities: 10% and 30%. The result is shown in Figure 11. It can be found that when the number of voxels reaches to 100, the result has a trend to converge. However, because the computational cost is quite large when the resolution is 100, there is no need to set such a high resolution if high accuracy is not the primary goal of simulation. A resolution around 50 can already give an indicative result. But if a high accuracy is required, 100 resolution is recommended.

Figure 11 The relationship between the resolution of the voxel and the axial Young's modulus

5.Extension

Shifted incomplete Cholesky factorization

To further accelerate the convergence of the PCG method, the shifted incomplete Cholesky factorization with threshold dropping can be used to modify the preconditioner. The default preconditioner in this code is an incomplete Cholesky factorization with zerofill. Nevertheless, the incomplete Cholesky with threshold dropping (ICT) can perform better convergence than zero-fill. A nonnegative scalar is used as a drop tolerance when performing ICT. Elements which are smaller in magnitude than a local drop tolerance are dropped from the resulting factor except for the diagonal element which is never dropped. However, Incomplete Cholesky factorizations of positive definite matrices do not always exist. Sometimes non-positive will be encountered when constructing an incomplete Cholesky preconditioner. To solve this problem, the 'diagcomp' option is used to construct a shifted incomplete Cholesky factorization. That is, instead of constructing L such that L*L' approximates A, incomplete Cholesky factorization with diagonal compensation constructs L such that approximates without explicitly forming M. As incomplete factorizations always exist for diagonally dominant matrices, alpha can be found to make M diagonally dominant. These can be done by substituting line 61 in Appendix A with the following lines: where 'opts' saves the options for 'ichol', the drop tolerance is 0.001 and the diagonal compensation is 0.01. To show the improvement of the shifted incomplete Cholesky factorization, a PCG method is performed on the Octet topology with a 10% relative density. The result is shown in Figure 12. It can be found that the result didn't converge to in 300 iterations when incomplete Cholesky preconditioner had zero-fill. When the threshold dropping is active and the diagonal compensation is 0.1, the result converged to in 250 iterations. Furthermore, when the diagonal compensation is 0.01, the result converged in 150 steps. But it is also found that if there is no diagonal compensation or the diagonal compensation is less than 0.001, the non-positive pivot error will be encountered. Therefore, a proper diagonal compensation should be selected in the shifted incomplete Cholesky factorization.

Two or more materials

Recently, Multi-material structures have received more and more attention from researchers, and have been successfully fabricated by some AM technologies [START_REF] Vaezi | Multiple material additive manufacturing -Part 1: a review[END_REF]. Cellular architectures fabricated with multi-materials can realize specific functions such as controllable thermal expansion [START_REF] Xu | Structurally efficient three-dimensional metamaterials with controllable thermal expansion[END_REF] and self-deformation [START_REF] Yu | Morpho: A self-deformable modular robot inspired by cellular structure[END_REF]. To simulate the material properties of multi-material cellular architectures, the homogenization code can be extended to assign more than one material properties. In this section, the extension to three materials is described. Firstly, the indicator matrix will not only have '0' and '1' elements to represent the void and the first material, but also have '2' and '3' elements to represent the second and the third material, respectively. Because the bulk material property is divided into two parts, and . If there is only one material, and are scalars. But if there are more than one materials, and will be vectors. Line 45 in Appendix A needs to be replaced by The vector and should have 3 entries, which are the Lame's first and second parameter for three materials, respectively. Also, the active degrees of freedom should also be changed to consider the elements with the second material and the third material.

The line 62 should be substituted by which means if the element has material 1, 2 or 3, its degrees of freedom will be activated in the global stiffness equation. Figure 13 shows an example of a three-material lattice structure model with 10% relative density. The Cubic-center topology can be disassembled to Cubic topology and X-shape topology. The voxels in the Cubic topology is assigned with material 1; the voxels in the X-shape topology is assigned with material 2; the connection between these two topologies is assigned with material 3. These three types of voxels are assembled together to get the multi-material Cubic-center topology.

Material 1 Material 2 Material 3 To get the homogenized constitutive matrix of this three-material structure, the material properties is assigned with the data in Table 3. The anisotropy can be controlled by adjusting the material properties. Calling the modified function as: The homogenized constitutive matrix for the multi-material lattice structure is obtained.

CH = homo3D(1,
To visualize the elastic properties along different orientations, the Young's modulus surface of the multi-materials lattice is shown in Figure 14(a). Figure 14(b) shows the Young's modulus surface of the single material lattice with and . The multi-materials lattice has an isotropic elastic property because the Young's modulus is the same along all the directions. Therefore, this homogenization code can help designers to design lattice structure with isotropic material properties.

Thermal conductivity

In thermal conductivity problems, the temperature is a scalar field. The homogenization equations are analogous to those of the elastic problem, which can be written as [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF]:

where T is the temperature field, is a virtual temperature field and is the conductivity tensor and. Unlike the elastic problem, only one material parameter is used to determine the conductivity of an isotropic material. Therefore, the input argument 'lambda' is a zero vector. The conductivity of different materials is given by the input argument 'mu'. To change it in the code, line 101 is substituted by: CMu = diag([1 1 1 0 0 0]); CLambda = zeros(6);

The element matrix 'keMu' contains the contributions in x, y and z directions in different rows and columns. Because only a scalar field is necessary, the contributions in x, y and z directions can be summed together to form an element matrix. The summation is done by adding the following line after line 15: Even though 'lambda' has no meaning in the thermal conductivity problem. 'lambda' is kept here so that the code can be modified as less as possible. Because every three rows and columns are saved in one row and column, the solution of the global stiffness equation should be modified. Line 59-66 is replaced by the following lines: L = ichol(K(4:3:end,4:3:end)); X = zeros(ndof,3); for i = 1:3 X(4:3:end,i) = pcg(K(4:3:end,4:3:end),F(3+i:3:end,i),1e-10,300,L,L'); end Also, element strain cases in line 75 should be modified as: The loops start from line 85 and line 85 should only run 3 times instead of 6 times. By modifying the code like this, the 3D conductivity problem with two materials can be solved. If three or more materials exist in the unit cell, the code can be further extended by the procedure similar to that in Section 5.2.

6.Conclusion

In this paper, a simple and efficient 3D homogenization code written in Matlab is provided to reduce the barrier of the implementation of numerical homogenization. Designers can use this code to obtain the homogenized constitutive matrix of 3D cellular materials or multi-material composites even without strong mathematical background. To further facilitate the implementation of this code, an algorithm is proposed to generate the voxel model of lattice structures. This model can be directly imported into the homogenization code. A wireframe script is needed in the generation of voxel models. Its format is defined so that users can customize the topology of the lattice structure.

By using the proposed method, the homogenized constitutive matrix and the Young's modulus surface of lattice structures with seven topologies are calculated. The results have shown that the lattice structures exhibit high anisotropy. The Cubic, Face-center, Tesseract, and Vintiles topologies have a higher stiffness along axial directions. On the contrary, the X-shape, Cubic-center and Octet topologies have a higher stiffness along diagonal directions.

The convergence of the PCG method in the homogenization code is accelerated by using Shifted incomplete Cholesky factorization. The code is then extended to multi-materials and thermal conductivity problems. It is found that the anisotropy of multi-material lattice structures can be controlled. With the proper combination of material properties, the Cubic center lattice structure is able to have isotropic elastic properties. Finally, future work can focus on the further improvement of the efficiency. For instance, the assembly-free method can be implemented to reduce the computational cost. Figure 12 The convergence plot for PCG method with different preconditioners ........... 17 List of Tables

List of Figures

Figure 1

 1 Figure 1 Wireframe of seven topologies. The next step is to use the information of the wireframe to generate the voxel model. The Matlab code provided in Appendix A named 'GenerateVoxel.m' computes the voxel model. The first argument (n) is the resolution of the voxel, which means the number of voxels along x, y and z-axis. The second argument (address) is the address of the wireframe file. The third argument (radius) is the radius of the lattice strut. There are two outputs. The first output 'voxel' is a 3D logical matrix, in which '0' means the voxel in this place does not have material; '1' means this voxel has material. The second output measures the relative density of the lattice structure.

Figure 2

 2 Figure 2 An example of the voxel model of a lattice structure with cubic topology.

Figure 3

 3 Figure 3 Two cases when calculating the minimum distance from voxel center to the line segment of lattice wireframe.

Figure 4

 4 Figure 4 Voxel models of seven topologies.

Figure 5

 5 Figure 5 Local node numbers and the natural coordinate of hexahedron element.

Figure 6

 6 Figure 6 An example of hexahedron mesh used to model the unit cell, (a) the ID of each element, (b) the non-periodic node number, (c) global node number of element 1 compared to local node number, (d) periodic node numbers.

Figure 7

 7 Figure 7 An example of periodic cellular materials with randomized unit cell

Figure 9 Figure 10

 910 Figure 9 Effective Young's modulus of lattice structures with 30% relative density.

  opts.type = 'ict'; opts.droptol = 1e-3; alpha = 0.01; opts.diagcomp = alpha; L = ichol(K(activedofs(4:end),activedofs(4:end)),opts);

Figure 12

 12 Figure 12 The convergence plot for PCG method with different preconditioners

  lambda = lambda(1)*(voxel==1) + lambda(2)*(voxel==2)... + lambda(3)*(voxel==3); mu = mu(1)*(voxel==1) + mu(2)*(voxel==2)... + mu(3)*(voxel==3);

  activedofs = edof((voxel==1 | voxel==2 | voxel==3),:);

Figure 13 A

 13 Figure 13 A three-material lattice cellular structure with Cubic-center topologies

Figure 14

 14 Figure 14 The Young's modulus surface of Cubic-center lattice with (a) multi-materials, (b) single material

keMu( 1 : 3 :

 13 end,1:3:end) = keMu(1:3:end,1:3:end)... + keMu(2:3:end,2:3:end) + keMu(3:3:end,3:3:end); To consider the conductivity of both materials in the cellular architecture, line 45 is changed to: lambda = lambda(1)*(voxel==0) + lambda(2)*(voxel==1); mu = mu(1)*(voxel==0) + mu(2)*(voxel==1);

X0_e( 4 : 3 :

 43 end,1:3) = keMu(4:3:end,4:3:end)... \[feMu(4:3:end,1),feMu(5:3:end,2),feMu(6:3:end,3)];

Figure 1 4 Figure 2 5 Figure 3 5 Figure 4 6 Figure 5 6 Figure 6 8 Figure 7 12 Figure 8 13 Figure 9

 1425354656687128139 Figure 1 Wireframe of seven topologies. ............................................................................ 4 Figure 2 An example of the voxel model of a lattice structure with cubic topology. ......... 5 Figure 3 Two cases when calculating the minimum distance from voxel center to the line segment of lattice wireframe. .............................................................................................. 5 Figure 4 Voxel models of seven topologies. ....................................................................... 6 Figure 5 Local node numbers and the natural coordinate of hexahedron element. ............ 6 Figure 6 An example of hexahedron mesh used to model the unit cell, (a) the ID of each element, (b) the non-periodic node number, (c) global node number of element 1 compared to local node number, (d) periodic node numbers. ............................................. 8 Figure 7 An example of periodic cellular materials with randomized unit cell ................ 12 Figure 8 Effective Young's modulus of lattice structures with 10% relative density. ..... 13 Figure 9 Effective Young's modulus of lattice structures with 30% relative density. ..... 13 Figure 10 Effective Young's modulus of lattice structures with 50% relative density. ... 14 Figure 11 The relationship between the resolution of the voxel and the axial Young's modulus ............................................................................................................................. 15 Figure12The convergence plot for PCG method with different preconditioners ........... 17 Figure 13 A three-material lattice cellular structure with Cubic-center topologies .......... 18 Figure 14 The Young's modulus surface of Cubic-center lattice with (a) multi-materials, (b) single material ............................................................................................................. 19

Figure 13 A

 13 Figure 1 Wireframe of seven topologies. ............................................................................ 4 Figure 2 An example of the voxel model of a lattice structure with cubic topology. ......... 5 Figure 3 Two cases when calculating the minimum distance from voxel center to the line segment of lattice wireframe. .............................................................................................. 5 Figure 4 Voxel models of seven topologies. ....................................................................... 6 Figure 5 Local node numbers and the natural coordinate of hexahedron element. ............ 6 Figure 6 An example of hexahedron mesh used to model the unit cell, (a) the ID of each element, (b) the non-periodic node number, (c) global node number of element 1 compared to local node number, (d) periodic node numbers. ............................................. 8 Figure 7 An example of periodic cellular materials with randomized unit cell ................ 12 Figure 8 Effective Young's modulus of lattice structures with 10% relative density. ..... 13 Figure 9 Effective Young's modulus of lattice structures with 30% relative density. ..... 13 Figure 10 Effective Young's modulus of lattice structures with 50% relative density. ... 14 Figure 11 The relationship between the resolution of the voxel and the axial Young's modulus ............................................................................................................................. 15 Figure12The convergence plot for PCG method with different preconditioners ........... 17 Figure 13 A three-material lattice cellular structure with Cubic-center topologies .......... 18 Figure 14 The Young's modulus surface of Cubic-center lattice with (a) multi-materials, (b) single material ............................................................................................................. 19

Table 1 The comparison between the PCG method and the direct method

 1

	PCG method

Table 2 Strut radius of lattice structures with different topologies

 2 Because the elastic properties of lattice structures exhibit high anisotropy, the homogenized Young's modulus varies in different orientations[START_REF] Xu | Design of lattice structures with controlled anisotropy[END_REF]. To illustrate the anisotropy of different types of lattice structures, the effective Young's modulus surfaces are used to display the Young's modulus in all the directions in a 3D space as shown in Figure8-10.

		Cubic	X-shape	Cubic	Face	Octet	Tesseract Vintiles
				center	center
	10%	0.113	0.074	0.063	0.058	0.045	0.054	0.065
	30%	0.208	0.134	0.114	0.109	0.087	0.100	0.128
	50%	0.287	0.187	0.157	0.150	0.120	0.138	0.177

Figure 8 Effective Young's modulus of lattice structures with 10% relative density.

			30% Relative Density
	Cubic	GPa	X-shape	GPa	Cubic center	GPa	Face center	GPa
		GPa			GPa		GPa
		Octet		Tesseract		Vintiles

Table 3 Material properties of three materials in the multi-material structure

 3

Table 1

 1 The comparison between the PCG method and the direct method ........................ 9 Table 2 Strut radius of lattice structures with different topologies ................................... 12 Table 3 Material properties of three materials in the multi-material structure ................. 18 Table A1 Definition and description of GRID………………………………….……….. 21 Table A2 Definition and description of STRUT………………………………………....21

Acknowledgment

This research work is supported by National Sciences and Engineering Research Council of Canada Discovery Grant RGPIN436055-2013.

Appendix A. Matlab files and the wireframe definition

All the Matlab files including the 'homo3D.m' and the 'GenerateVoxel.m' and other supporting files such as the wireframe files can be downloaded by the link: https://github.com/GuoyingDong/homogenization or https://www.mathworks.com/matlabcentral/fileexchange/67457-3d-homogenization-ofcellular-materials