Metal oxide sunscreen

Note: the cost of running the entire notebook is larger than 1 FlexCredits.

Sunscreen technology, developed in the 1940s to mitigate sunburn, has advanced significantly with the integration of metal oxide particles, primarily zinc oxide and titanium dioxide. These formulations provide robust protection against UV rays, ensuring effective defense against varying intensities of sun exposure. A common misconception is that the protective action mainly arises from the reflection and scattering of these active metal oxide particles, rather than their absorption.

In this notebook, we will reproduce the results presented in the paper titled: Cole, C., Shyr, T., & Ou‐Yang, H. (2015). Metal oxide sunscreens protect skin by absorption, not by reflection or scattering https://doi.org/10.1111/phpp.12214, where the authors confirmed, through absorption and reflectance measurements, that the main mechanism underlying UV block in metal oxide particle-based sunscreen relies on absorption of the particles, rather than reflection or scattering.

The simulation model consists of an agglomerate of TiO₂ nanoparticles excited by a plane wave at normal incidence. Flux monitors are positioned on the top and bottom of the simulation domain to compute transmitted and forward scattered fields (which would penetrate the skin) and reflected and back-scattered fields (which do not penetrate the skin due to the nanoparticles action). The material properties of the TiO₂ particles are obtained directly from refractiveindex.info, and fitted using Tidy3D’s dispersion fitting tool.

For more examples involving scattering of nanoparticles and structured surfaces, please refer to our example library, where you can find interesting case studies such as scattering cross-section calculation of a dielectric sphere, the scattering of a plasmonic nanoparticle, the gradient metasurface reflector, and the dielectric metasurface absorber.

Simulation Setup

First, we will import the necessary libraries.

import tidy3d as td
import numpy as np
import matplotlib.pyplot as plt
from tidy3d import web

# defining a random seed for reproducibility
np.random.seed(12)

Retrieving refractive index data and creating the medium.

We will retrieve material property data from refractiveindex.info using the from_url method of the `FastDispersionFitter <https://docs.flexcompute.com/projects/tidy3d/en/v2.5.2/_autosummary/tidy3d.plugins.dispersion.FastDispersionFitter.html>`__ class, and fit a pole residue medium.

For more information about fitting dispersive material, please refer to this example notebook.

from tidy3d.plugins.dispersion import FastDispersionFitter

url = 'https://refractiveindex.info/tmp/database/data-nk/main/TiO2/Zhukovsky.txt'

# creating the fitter object loading that from the url
fitter = FastDispersionFitter.from_url(url,delimiter= '\t')

# fitting the data
medium, error = fitter.fit(max_num_poles=5)

13:19:38 -03 WARNING: Unable to fit with weighted RMS error under               
             'tolerance_rms' of 1e-05                                           

Although the RMS error cannot be below \(10^{-5}\), we can see that the material properties are well described by the fit.

# plotting for comparison
fig, ax = plt.subplots()

# original dataset
freq = td.C_0 / fitter.wvl_um * 10**-12
ax.plot(freq, fitter.n_data, 'o', alpha=0.2, color="black", label = 'n')
ax.plot(freq, fitter.k_data, 'o', alpha=0.2, color="black", label = 'k')


medium.plot(np.linspace(td.C_0 / 0.215, td.C_0 / 1.5, 1000), ax=ax)
ax.lines[-1].set_ls("--")
ax.lines[-1].set_label('Fitted n')
ax.lines[-2].set_ls("--")
ax.lines[-2].set_label('Fitted k')

Define parameters for building the simulation

We are now defining the parameters that will be used to build the simulation.

# for simplicity, the nanoparticles are modeled as spheres
size = 0.025

# volume to place the nanoparticles
lx = 1.2
ly = 1.2
lz = 10 * size

# target spectrum
wl1 = 0.22
wl2 = 0.800

central_wl = (wl2 - wl1) / 2 + wl1
freq1 = td.C_0 / wl1
freq2 = td.C_0 / wl2
freq0 = (freq1 - freq2) / 2 + freq2
fwidth = (freq1 - freq2) / 2

# frequencies to be analyzed in the monitors
freqs = np.linspace(freq0 - fwidth, freq0 + fwidth, 1000)

# size of the simulation domain
Lz = 2 * wl2 + lz
Lx = lx
Ly = ly

# souce and monitor position
monitor_z = (Lz / 2 - lz / 2) / 3
source_z = -Lz / 2 + 2 * (Lz / 2 - lz / 2) / 3

structures_center = (0, 0, 0)

# simulation runtime
run_time = 5e-12

# boundary specifications. We will define PML at the z-boundaries and periodic conditions on the remaining ones."
boundary_spec = td.BoundarySpec(
    x=td.Boundary.periodic(), y=td.Boundary.periodic(), z=td.Boundary.pml()
)

Now we will create an auxiliary function to randomly place the nanoparticles as a function of the concentration. The number of particles is calculated with parameters given in the paper.

def get_particles(proportion=0.2):
    # concentration of petrolatum matrix
    petrolatum_concentration = 1.3e-8  # mg/um^3
    active_proportion = proportion
    particle_volume = (4 / 3) * np.pi * size**3
    Ti02_density = 4.23e-9  # mg/um^3
    particle_weight = Ti02_density * particle_volume

    # calculate the total mass of Ti02 nanoparticles
    volume = lx * ly * lz
    Ti02_mass = active_proportion * volume * petrolatum_concentration  # mg

    # defining the number of particles
    N_particles = int(Ti02_mass / particle_weight)

    # create arrays of random positions for each particle
    positions_x = np.random.uniform(-lx / 2, lx / 2, N_particles)
    position_y = np.random.uniform(-ly / 2, ly / 2, N_particles)
    positions_z = np.random.uniform(
        structures_center[2] - lz / 2, structures_center[2] + lz / 2, N_particles
    )

    # create a first geometry object and sum the others to form the agglomerate
    particles = td.Sphere(
        center=(positions_x[0], position_y[0], positions_z[0]), radius=size
    )

    for pos in zip(positions_x[1:], position_y[1:], positions_z[1:]):
        particles += td.Sphere(center=pos, radius=size)

    particles = td.Structure(geometry=particles, medium=medium)

    return [particles]

Define source and monitors:

The source spectrum is chosen to cover the most intense part of the solar emission and the UV region, which is the main concern for skin health.

Since the size of the nanoparticles is much smaller than the target wavelengths, we will add a mesh override structure to properly resolve the nanoparticles agglomerate.

We will simulate for two active particle concentration, 0.1% and 0.2%.

source = td.PlaneWave(
    center=(0, 0, source_z),
    size=(td.inf, td.inf, 0),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    direction="+",
)

monitor_bottom = td.FluxMonitor(
        center=(0, 0, -Lz / 2 + monitor_z),
        size=(td.inf, td.inf, 0),
        freqs=freqs,
        name="z-",
    )

monitor_top = td.FluxMonitor(
        center=(0, 0, Lz / 2 - monitor_z),
        size=(td.inf, td.inf, 0),
        freqs=freqs,
        name="z",
    )



mesh_override = td.MeshOverrideStructure(
    geometry=td.Box(center=structures_center, size=(lx, ly, lz + 2 * size)),
    dl=(size / 6,) * 3,
)
grid_spec = td.GridSpec.auto(min_steps_per_wvl=20, override_structures=[mesh_override])

sim = td.Simulation(
    size=(Lx, Ly, Lz),
    structures=[],
    sources=[source],
    monitors=[monitor_bottom,monitor_top],
    run_time=run_time,
    boundary_spec=boundary_spec,
    grid_spec=grid_spec,
)

Simulations = {
    "proportion 0.1%": sim.updated_copy(structures=get_particles(0.1)),
    "proportion 0.2%": sim.updated_copy(structures=get_particles(0.2)),
}

# visualizing
Simulations["proportion 0.1%"].plot_3d()

We will run both simulations in parallel using the `web.Batch <https://docs.flexcompute.com/projects/tidy3d/en/v1.5.0/_autosummary/tidy3d.web.container.Batch.html>`__ function.

batch = web.Batch(simulations=Simulations, verbose=True)
results = batch.run(path_dir="batch_simulations")

13:19:52 -03 Started working on Batch containing 2 tasks.

13:19:55 -03 Maximum FlexCredit cost: 16.169 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.

13:26:14 -03 Batch complete.

We will organize the data in a pandas DataFrame object.

The transmitted flux is defined as the flux in the upper monitor, the reflected flux as the flux in the bottom monitor, and the absorbed flux will be defined as \(A = 1 - (T + R)\)

import pandas as pd

df = pd.DataFrame(
    index=td.C_0 / freqs,
    columns=[
        "Transmitted 0.1%",
        "Reflected 0.1%",
        "Absorbed 0.1%",
        "Transmitted 0.2%",
        "Reflected 0.2%",
        "Absorbed 0.2%",
    ],
)

for sim in ["proportion 0.1%", "proportion 0.2%"]:
    sim_data = results[sim]

    T = sim_data["z"].flux
    R = abs(sim_data["z-"].flux)
    A = 1 - (T + R)

    df["Transmitted %s" % sim.split()[-1]] = T
    df["Reflected %s" % sim.split()[-1]] = R
    df["Absorbed %s" % sim.split()[-1]] = A

df.head(10)

	Transmitted 0.1%	Reflected 0.1%	Absorbed 0.1%	Transmitted 0.2%	Reflected 0.2%	Absorbed 0.2%
0.800000	0.976410	0.032138	-0.008548	0.929303	0.010096	0.060601
0.797894	0.976502	0.031542	-0.008044	0.929808	0.009658	0.060535
0.795800	0.976489	0.030921	-0.007410	0.930475	0.009248	0.060277
0.793716	0.976414	0.030284	-0.006697	0.930902	0.008905	0.060194
0.791643	0.976394	0.029661	-0.006055	0.931077	0.008639	0.060284
0.789581	0.976473	0.029066	-0.005539	0.931135	0.008458	0.060407
0.787530	0.976561	0.028478	-0.005039	0.930816	0.008351	0.060833
0.785490	0.976556	0.027873	-0.004429	0.930435	0.008318	0.061247
0.783460	0.976475	0.027251	-0.003726	0.929848	0.008350	0.061803
0.781440	0.976435	0.026638	-0.003073	0.929037	0.008442	0.062521

Plotting the results, we can see the transmittance follows a similar behavior as reported for the TiO₂ < 25 nm nanoparticles in Figure 2 of the paper.

Looking at the proportion of reflected light, we can conclude that the main mechanism of UV blocking is indeed the absorption of the nanoparticles.

It is interesting to note the complete blockage in the UVB region (280 - 315 nm) and significant attenuation in the UVA region (315 - 400 nm), especially for the 0.2% proportion, primarily due to absorption.

fig, ax = plt.subplots()
ax = df["Transmitted 0.1%"].plot()
ax = df["Reflected 0.1%"].plot(ax=ax)
ax = df["Absorbed 0.1%"].plot(ax=ax)
ax.set_xlabel('wavelength ($\mu$m)')
ax.legend()

fig, ax = plt.subplots()
ax = df["Transmitted 0.2%"].plot()
ax = df["Reflected 0.2%"].plot(ax=ax)
ax = df["Absorbed 0.2%"].plot(ax=ax)
ax.set_xlabel('wavelength ($\mu$m)')
ax.legend()

plt.show()

print(
    "Proportion of scattered light at 0.1%s: %.2f" % ("%", df["Reflected 0.1%"].mean())
)
print(
    "Proportion of scattered light at 0.2%s: %.2f" % ("%", df["Reflected 0.2%"].mean())
)

Proportion of scattered light at 0.1%: 0.05
Proportion of scattered light at 0.2%: 0.08