Chapter 4: Functional Ceramics

4.1 Overview of Functional Ceramics

Functional ceramics are material groups that utilize dielectric, piezoelectric, ionic conduction, luminescence, and magnetic properties. They are widely applied as capacitors, sensors, solid electrolytes, LED phosphors, and magnetic devices. In this chapter, we will learn the physical mechanisms of each function and practice property calculations with Python.

Classification of Functional Ceramics

4.2 Dielectric Ceramics

4.2.1 Mechanism of Dielectric Polarization

When an electric field is applied to a dielectric, the following four types of polarization occur:

1. Electronic polarization: Displacement of electron cloud (all materials, ~10¹⁵ Hz)
2. Ionic polarization: Relative displacement of ions (infrared range, ~10¹³ Hz)
3. Orientation polarization: Rotation of permanent dipoles (microwave range, ~10⁹ Hz)
4. Interfacial polarization: Charge accumulation at interfaces (low frequency, ~10³ Hz)

The permittivity \( \epsilon_r \) is defined as the ratio to the vacuum permittivity \( \epsilon_0 = 8.854 \times 10^{-12} \) F/m:

\[ \epsilon_r = \frac{\epsilon}{\epsilon_0} \]

4.2.2 Ferroelectricity of BaTiO₃

BaTiO₃ undergoes a phase transition to tetragonal structure below the Curie temperature (T_c ≈ 120°C) and possesses spontaneous polarization. This ferroelectricity results in extremely high permittivity values exceeding 10,000.

Python Implementation: Permittivity Calculation using Clausius-Mossotti Equation

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 1: Clausius-Mossotti Permittivity Calculation
# ===================================

import numpy as np
import matplotlib.pyplot as plt

def clausius_mossotti_permittivity(alpha, N, epsilon_0=8.854e-12):
    """
    Calculate permittivity using the Clausius-Mossotti equation

    Parameters:
    -----------
    alpha : float
        Polarizability [C·m²/V]
    N : float
        Number of atoms per unit volume [m^-3]
    epsilon_0 : float
        Vacuum permittivity [F/m]

    Returns:
    --------
    epsilon_r : float
        Relative permittivity (dimensionless)
    """
    # Clausius-Mossotti equation: (εᵣ - 1)/(εᵣ + 2) = N·α/(3ε₀)
    # Rearranged: εᵣ = (1 + 2·N·α/(3ε₀)) / (1 - N·α/(3ε₀))

    numerator = 1 + 2 * N * alpha / (3 * epsilon_0)
    denominator = 1 - N * alpha / (3 * epsilon_0)

    epsilon_r = numerator / denominator

    return epsilon_r


def curie_weiss_law(T, T_c=393, C=1.5e5):
    """
    Temperature dependence of permittivity using Curie-Weiss law

    Parameters:
    -----------
    T : float or array
        Temperature [K]
    T_c : float
        Curie temperature [K] (BaTiO₃ ≈ 120°C = 393K)
    C : float
        Curie constant [K]

    Returns:
    --------
    epsilon_r : float or array
        Relative permittivity
    """
    # Curie-Weiss law: εᵣ = C / (T - Tᶜ)
    epsilon_r = C / np.abs(T - T_c)

    # Limit maximum value to avoid divergence
    epsilon_r = np.minimum(epsilon_r, 20000)

    return epsilon_r


def plot_dielectric_properties():
    """
    Visualize dielectric properties
    """
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Polarizability vs permittivity relationship
    ax1 = axes[0]

    # Typical material polarizability range
    alpha_range = np.linspace(1e-40, 5e-40, 100)  # C·m²/V
    N = 5e28  # m^-3 (typical ion density)

    epsilon_r_values = [clausius_mossotti_permittivity(alpha, N) for alpha in alpha_range]

    ax1.plot(alpha_range * 1e40, epsilon_r_values, linewidth=2, color='navy')

    # Plot representative materials
    materials = {
        'SiO₂': {'alpha': 3.5e-41, 'epsilon_r': 3.8},
        'Al₂O₃': {'alpha': 1.0e-40, 'epsilon_r': 9.0},
        'TiO₂': {'alpha': 2.5e-40, 'epsilon_r': 80},
        'BaTiO₃': {'alpha': 4.5e-40, 'epsilon_r': 1200}
    }

    for name, props in materials.items():
        ax1.scatter(props['alpha'] * 1e40, props['epsilon_r'],
                   s=150, edgecolors='black', linewidth=2, zorder=5)
        ax1.annotate(name, (props['alpha'] * 1e40, props['epsilon_r']),
                    xytext=(5, 5), textcoords='offset points', fontsize=10)

    ax1.set_xlabel('Polarizability α (10^-40 C·m²/V)', fontsize=12)
    ax1.set_ylabel('Relative Permittivity εᵣ', fontsize=12)
    ax1.set_title('Clausius-Mossotti Relation', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)

    # Right plot: BaTiO₃ Curie-Weiss behavior
    ax2 = axes[1]

    T_range = np.linspace(293, 493, 500)  # 20°C ~ 220°C
    T_c = 393  # 120°C

    epsilon_r_bto = curie_weiss_law(T_range, T_c, C=1.5e5)

    ax2.plot(T_range - 273.15, epsilon_r_bto, linewidth=2, color='crimson')
    ax2.axvline(x=T_c - 273.15, color='blue', linestyle='--', linewidth=2,
                label=f'Curie temperature = {T_c-273.15:.0f}°C')

    ax2.set_xlabel('Temperature (°C)', fontsize=12)
    ax2.set_ylabel('Relative Permittivity εᵣ', fontsize=12)
    ax2.set_title('BaTiO₃ Curie-Weiss Behavior', fontsize=14, fontweight='bold')
    ax2.set_ylim(0, 20000)
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('dielectric_properties.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Numerical examples
    print("=== Clausius-Mossotti Calculation ===")
    for name, props in materials.items():
        calculated_eps = clausius_mossotti_permittivity(props['alpha'], N)
        print(f"{name:10s}: α = {props['alpha']*1e40:.2f}×10^-40, εᵣ = {calculated_eps:.1f} (exp: {props['epsilon_r']})")

    print("\n=== BaTiO₃ Temperature Dependence ===")
    for T_celsius in [80, 100, 120, 140, 160]:
        T = T_celsius + 273.15
        eps = curie_weiss_law(T, T_c=393, C=1.5e5)
        print(f"T = {T_celsius:3d}°C → εᵣ = {eps:6.0f}")

# Execute
plot_dielectric_properties()

4.2.3 MLCC (Multilayer Ceramic Capacitor) Design

MLCCs achieve high capacitance through a structure of alternating dielectric layers and internal electrodes. The capacitance C is:

\[ C = \frac{\epsilon_0 \epsilon_r A}{d} \times n \]

where \( A \) is the electrode area, \( d \) is the dielectric layer thickness, and \( n \) is the number of layers.

Python Implementation: MLCC Capacitance Design Calculation

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 2: MLCC Capacitance Design
# ===================================

import numpy as np
import matplotlib.pyplot as plt

def mlcc_capacitance(epsilon_r, A, d, n, epsilon_0=8.854e-12):
    """
    Calculate capacitance of MLCC (Multilayer Ceramic Capacitor)

    Parameters:
    -----------
    epsilon_r : float
        Relative permittivity
    A : float
        Electrode area [m²]
    d : float
        Dielectric layer thickness [m]
    n : int
        Number of layers
    epsilon_0 : float
        Vacuum permittivity [F/m]

    Returns:
    --------
    C : float
        Capacitance [F]
    """
    C = (epsilon_0 * epsilon_r * A / d) * n
    return C


def design_mlcc_optimization():
    """
    MLCC design optimization simulation
    """
    # Design parameters
    epsilon_r = 3000  # X7R material (BaTiO₃-based)
    A = (1e-3)**2  # 1mm × 1mm
    d_range = np.linspace(1e-6, 50e-6, 100)  # 1~50 μm
    n = 100  # 100 layers

    # Capacitance calculation
    C_values = [mlcc_capacitance(epsilon_r, A, d, n) for d in d_range]

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Layer thickness vs capacitance
    ax1 = axes[0]
    ax1.plot(d_range * 1e6, np.array(C_values) * 1e6, linewidth=2, color='navy')
    ax1.axhline(y=1, color='red', linestyle='--', alpha=0.5, label='Target: 1 μF')

    ax1.set_xlabel('Dielectric Layer Thickness (μm)', fontsize=12)
    ax1.set_ylabel('Capacitance (μF)', fontsize=12)
    ax1.set_title('MLCC Capacitance vs Layer Thickness', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # Right plot: Effect of number of layers
    ax2 = axes[1]
    d_fixed = 5e-6  # 5 μm
    n_range = np.arange(10, 501, 10)
    C_vs_n = [mlcc_capacitance(epsilon_r, A, d_fixed, n) for n in n_range]

    ax2.plot(n_range, np.array(C_vs_n) * 1e6, linewidth=2, color='green')

    ax2.set_xlabel('Number of Layers n', fontsize=12)
    ax2.set_ylabel('Capacitance (μF)', fontsize=12)
    ax2.set_title(f'MLCC Capacitance vs Layer Count (d = {d_fixed*1e6} μm)', fontsize=14, fontweight='bold')
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('mlcc_design.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Design examples
    print("=== MLCC Design Examples ===")
    designs = [
        {'name': 'Low-C (1 nF)', 'd': 20e-6, 'n': 10},
        {'name': 'Mid-C (100 nF)', 'd': 10e-6, 'n': 50},
        {'name': 'High-C (1 μF)', 'd': 5e-6, 'n': 200},
        {'name': 'Ultra-C (10 μF)', 'd': 2e-6, 'n': 500}
    ]

    for design in designs:
        C = mlcc_capacitance(epsilon_r, A, design['d'], design['n'])
        volume = A * design['d'] * design['n']
        density = C / volume  # F/m³
        print(f"{design['name']:20s}: d={design['d']*1e6:4.1f} μm, n={design['n']:3d} → C={C*1e6:8.2f} μF, density={density*1e-9:.2f} μF/mm³")

# Execute
design_mlcc_optimization()

4.3 Piezoelectric Ceramics

4.3.1 Principle of Piezoelectric Effect

The piezoelectric effect is the phenomenon where charge is generated by mechanical stress (direct piezoelectric effect) and strain is induced by applying an electric field (converse piezoelectric effect). Crystal structure asymmetry is necessary, and 20 of the 32 crystal point groups exhibit piezoelectricity.

The piezoelectric constant \( d_{33} \) represents the strain (or charge) when the electric field direction and stress direction are parallel:

\[ S_3 = d_{33} E_3 \quad \text{(Converse piezoelectric effect)} \] \[ D_3 = d_{33} T_3 \quad \text{(Direct piezoelectric effect)} \]

where \( S \) is strain, \( E \) is electric field, \( D \) is electric displacement, and \( T \) is stress.

4.3.2 Properties of PZT (Lead Zirconate Titanate)

PZT (Pb(Zr_xTi_1-x)O₃) is the most widely used piezoelectric material. Piezoelectric properties are maximized at the MPB (Morphotropic Phase Boundary, x ≈ 0.52) composition.

Python Implementation: Calculation of Piezoelectric Constant and Coupling Coefficient

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 3: Piezoelectric Constant and Coupling Coefficient Calculation
# ===================================

import numpy as np
import matplotlib.pyplot as plt

def piezoelectric_strain(d_33, E_3):
    """
    Calculate strain from converse piezoelectric effect

    Parameters:
    -----------
    d_33 : float
        Piezoelectric constant [m/V] or [C/N]
    E_3 : float
        Electric field [V/m]

    Returns:
    --------
    S_3 : float
        Strain (dimensionless)
    """
    S_3 = d_33 * E_3
    return S_3


def electromechanical_coupling_factor(d_33, epsilon_33, s_33):
    """
    Calculate electromechanical coupling coefficient

    Parameters:
    -----------
    d_33 : float
        Piezoelectric constant [C/N]
    epsilon_33 : float
        Permittivity [F/m]
    s_33 : float
        Elastic compliance [m²/N]

    Returns:
    --------
    k_33 : float
        Coupling coefficient (dimensionless, 0~1)
    """
    k_33 = d_33 / np.sqrt(epsilon_33 * s_33)
    return k_33


def plot_piezoelectric_response():
    """
    Visualize piezoelectric response
    """
    # Physical properties of PZT-5H
    d_33 = 593e-12  # C/N
    epsilon_0 = 8.854e-12
    epsilon_r = 3400
    epsilon_33 = epsilon_0 * epsilon_r

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Electric field-strain relationship (converse piezoelectric effect)
    ax1 = axes[0]
    E_range = np.linspace(-2e6, 2e6, 500)  # -2 ~ +2 MV/m
    strain = piezoelectric_strain(d_33, E_range)

    ax1.plot(E_range / 1e6, strain * 1e6, linewidth=2, color='navy')
    ax1.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
    ax1.axvline(x=0, color='black', linestyle='-', linewidth=0.5)

    ax1.set_xlabel('Electric Field E₃ (MV/m)', fontsize=12)
    ax1.set_ylabel('Strain S₃ (μstrain)', fontsize=12)
    ax1.set_title('Converse Piezoelectric Effect (PZT-5H)', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)

    # Right plot: Material comparison (piezoelectric constant and coupling coefficient)
    ax2 = axes[1]

    materials_piezo = {
        'BaTiO₃': {'d33': 190e-12, 'epsilon_r': 1700, 's33': 8.3e-12},
        'PZT-4': {'d33': 289e-12, 'epsilon_r': 1300, 's33': 12.3e-12},
        'PZT-5H': {'d33': 593e-12, 'epsilon_r': 3400, 's33': 16.5e-12},
        'PVDF': {'d33': 33e-12, 'epsilon_r': 12, 's33': 400e-12},
        'AlN': {'d33': 5.5e-12, 'epsilon_r': 9, 's33': 3.0e-12}
    }

    names = list(materials_piezo.keys())
    d33_values = [props['d33'] * 1e12 for props in materials_piezo.values()]
    k33_values = []

    for props in materials_piezo.values():
        epsilon = epsilon_0 * props['epsilon_r']
        k33 = electromechanical_coupling_factor(props['d33'], epsilon, props['s33'])
        k33_values.append(k33)

    x = np.arange(len(names))
    width = 0.35

    ax2.bar(x - width/2, d33_values, width, label='d₃₃ (pC/N)', color='skyblue', edgecolor='black')
    ax2_twin = ax2.twinx()
    ax2_twin.bar(x + width/2, k33_values, width, label='k₃₃', color='salmon', edgecolor='black')

    ax2.set_xlabel('Material', fontsize=12)
    ax2.set_ylabel('Piezoelectric Constant d₃₃ (pC/N)', fontsize=12, color='skyblue')
    ax2_twin.set_ylabel('Coupling Coefficient k₃₃', fontsize=12, color='salmon')
    ax2.set_title('Piezoelectric Material Comparison', fontsize=14, fontweight='bold')
    ax2.set_xticks(x)
    ax2.set_xticklabels(names, rotation=15)
    ax2.tick_params(axis='y', labelcolor='skyblue')
    ax2_twin.tick_params(axis='y', labelcolor='salmon')
    ax2.grid(True, alpha=0.3, axis='y')

    plt.tight_layout()
    plt.savefig('piezoelectric_response.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Numerical examples
    print("=== Piezoelectric Response Calculation (PZT-5H) ===")
    for E_MV in [0.5, 1.0, 2.0]:
        E = E_MV * 1e6
        S = piezoelectric_strain(d_33, E)
        displacement = S * 1e-3  # For 1 mm thickness
        print(f"E = {E_MV:.1f} MV/m → Strain = {S*1e6:.1f} μstrain, Displacement (1mm) = {displacement*1e9:.2f} nm")

    print("\n=== Electromechanical Coupling Coefficients ===")
    for name, props in materials_piezo.items():
        epsilon = epsilon_0 * props['epsilon_r']
        k33 = electromechanical_coupling_factor(props['d33'], epsilon, props['s33'])
        print(f"{name:10s}: d₃₃ = {props['d33']*1e12:5.1f} pC/N, k₃₃ = {k33:.3f}")

# Execute
plot_piezoelectric_response()

Python Implementation: PZT Phase Diagram and MPB Composition

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 4: PZT Phase Diagram and MPB Composition Visualization
# ===================================

import numpy as np
import matplotlib.pyplot as plt

def pzt_phase_diagram():
    """
    Visualize PZT phase diagram and MPB (Morphotropic Phase Boundary)
    """
    # Composition range (Zr content)
    x_Zr = np.linspace(0, 1, 500)

    # Composition dependence of Curie temperature (simplified model)
    T_c = 490 - 200 * x_Zr + 150 * x_Zr**2  # °C

    # Composition dependence of piezoelectric constant (maximum near MPB)
    # MPB composition: x ≈ 0.52
    x_MPB = 0.52
    d_33 = 200 + 400 * np.exp(-((x_Zr - x_MPB) / 0.1)**2)

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Phase diagram
    ax1 = axes[0]

    # Phase boundary definition (simplified)
    x_rhomb = x_Zr[x_Zr > x_MPB]
    x_tetra = x_Zr[x_Zr <= x_MPB]
    T_c_rhomb = T_c[x_Zr > x_MPB]
    T_c_tetra = T_c[x_Zr <= x_MPB]

    # Fill regions
    ax1.fill_between(x_rhomb, 0, T_c_rhomb, alpha=0.3, color='blue', label='Rhombohedral (FE_R)')
    ax1.fill_between(x_tetra, 0, T_c_tetra, alpha=0.3, color='red', label='Tetragonal (FE_T)')
    ax1.fill_between(x_Zr, T_c, 500, alpha=0.2, color='gray', label='Cubic (Paraelectric)')

    # MPB line
    ax1.axvline(x=x_MPB, color='green', linestyle='--', linewidth=2, label=f'MPB (x = {x_MPB})')
    ax1.plot(x_Zr, T_c, 'k-', linewidth=2, label='Curie Temperature')

    ax1.set_xlabel('Zr Content x in Pb(Zr_xTi_{1-x})O₃', fontsize=12)
    ax1.set_ylabel('Temperature (°C)', fontsize=12)
    ax1.set_title('PZT Phase Diagram', fontsize=14, fontweight='bold')
    ax1.set_xlim(0, 1)
    ax1.set_ylim(0, 500)
    ax1.legend(loc='upper right')
    ax1.grid(True, alpha=0.3)

    # Right plot: Composition dependence of piezoelectric constant
    ax2 = axes[1]
    ax2.plot(x_Zr, d_33, linewidth=2, color='navy')
    ax2.axvline(x=x_MPB, color='green', linestyle='--', linewidth=2, label=f'MPB (x = {x_MPB})')
    ax2.scatter([x_MPB], [np.max(d_33)], s=200, c='red', marker='*', edgecolors='black',
                linewidth=2, zorder=5, label=f'd₃₃ max = {np.max(d_33):.0f} pC/N')

    ax2.set_xlabel('Zr Content x in Pb(Zr_xTi_{1-x})O₃', fontsize=12)
    ax2.set_ylabel('Piezoelectric Constant d₃₃ (pC/N)', fontsize=12)
    ax2.set_title('PZT Piezoelectric Constant vs Composition', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('pzt_phase_diagram.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Numerical examples
    print("=== PZT Composition Effects ===")
    compositions = [0.3, 0.4, 0.52, 0.6, 0.7]
    for x in compositions:
        idx = np.argmin(np.abs(x_Zr - x))
        T_c_val = T_c[idx]
        d_33_val = d_33[idx]
        phase = 'Tetragonal' if x < x_MPB else 'Rhombohedral'
        print(f"x = {x:.2f} ({phase:15s}): T_c = {T_c_val:5.1f}°C, d₃₃ = {d_33_val:5.0f} pC/N")

# Execute
pzt_phase_diagram()

4.4 Ionic Conductor Ceramics

4.4.1 Mechanism of Ionic Conduction

In ionic conductors, ions (O²⁻, Li⁺, Na⁺, etc.) move through the crystal lattice. The conductivity \( \sigma \) is expressed by the Arrhenius equation:

\[ \sigma = \sigma_0 \exp\left(-\frac{E_a}{k_B T}\right) \]

where \( \sigma_0 \) is the pre-exponential factor, \( E_a \) is the activation energy, \( k_B \) is the Boltzmann constant, and \( T \) is the temperature.

4.4.2 YSZ (Yttria-Stabilized Zirconia)

YSZ ((Y₂O₃)_x(ZrO₂)_1-x) is used as the electrolyte in SOFCs (Solid Oxide Fuel Cells). Y³⁺ doping creates oxygen vacancies, allowing O²⁻ ion conduction.

Python Implementation: Arrhenius Conductivity and Activation Energy

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 5: Arrhenius Conductivity Calculation
# ===================================

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def ionic_conductivity(T, sigma_0, E_a, k_B=8.617e-5):
    """
    Calculate ionic conductivity using Arrhenius equation

    Parameters:
    -----------
    T : float or array
        Temperature [K]
    sigma_0 : float
        Pre-exponential factor [S/cm]
    E_a : float
        Activation energy [eV]
    k_B : float
        Boltzmann constant [eV/K]

    Returns:
    --------
    sigma : float or array
        Ionic conductivity [S/cm]
    """
    sigma = sigma_0 * np.exp(-E_a / (k_B * T))
    return sigma


def fit_arrhenius_data(T_data, sigma_data):
    """
    Fit Arrhenius parameters from measured data

    Parameters:
    -----------
    T_data : array
        Temperature data [K]
    sigma_data : array
        Conductivity data [S/cm]

    Returns:
    --------
    sigma_0 : float
        Pre-exponential factor [S/cm]
    E_a : float
        Activation energy [eV]
    """
    # Logarithmic transformation: ln(σ) = ln(σ₀) - E_a/(k_B·T)
    # y = ln(σ), x = 1/T, slope = -E_a/k_B, intercept = ln(σ₀)

    x = 1 / T_data
    y = np.log(sigma_data)

    # Linear regression
    coeffs = np.polyfit(x, y, 1)
    slope = coeffs[0]
    intercept = coeffs[1]

    E_a = -slope * 8.617e-5  # eV
    sigma_0 = np.exp(intercept)

    return sigma_0, E_a


def plot_ionic_conductivity():
    """
    Arrhenius plot and material comparison of ionic conductivity
    """
    # Parameters for each material
    materials_ionic = {
        'YSZ (8YSZ)': {'sigma_0': 2.5e4, 'E_a': 1.0, 'color': 'navy'},
        'GDC': {'sigma_0': 1.0e4, 'E_a': 0.7, 'color': 'green'},
        'LLZO': {'sigma_0': 1.0e2, 'E_a': 0.3, 'color': 'red'},
        'β-Al₂O₃': {'sigma_0': 5.0e2, 'E_a': 0.16, 'color': 'purple'}
    }

    # Temperature range
    T_range = np.linspace(300, 1200, 500)  # K

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Arrhenius plot
    ax1 = axes[0]

    for name, params in materials_ionic.items():
        sigma = ionic_conductivity(T_range, params['sigma_0'], params['E_a'])
        ax1.semilogy(1e4 / T_range, sigma, linewidth=2, label=name, color=params['color'])

    ax1.set_xlabel('10^4 / T (K^-1)', fontsize=12)
    ax1.set_ylabel('Ionic Conductivity σ (S/cm)', fontsize=12)
    ax1.set_title('Arrhenius Plot for Ionic Conductors', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3, which='both')

    # Right plot: Conductivity comparison at actual temperatures
    ax2 = axes[1]

    temperatures_celsius = [400, 600, 800, 1000]
    x_pos = np.arange(len(temperatures_celsius))
    width = 0.2

    for i, (name, params) in enumerate(materials_ionic.items()):
        sigma_values = [ionic_conductivity(T + 273.15, params['sigma_0'], params['E_a'])
                       for T in temperatures_celsius]
        ax2.bar(x_pos + i * width, sigma_values, width, label=name, color=params['color'], edgecolor='black')

    ax2.set_xlabel('Temperature (°C)', fontsize=12)
    ax2.set_ylabel('Ionic Conductivity σ (S/cm)', fontsize=12)
    ax2.set_title('Conductivity Comparison at Different Temperatures', fontsize=14, fontweight='bold')
    ax2.set_xticks(x_pos + width * 1.5)
    ax2.set_xticklabels(temperatures_celsius)
    ax2.set_yscale('log')
    ax2.legend()
    ax2.grid(True, alpha=0.3, which='both')

    plt.tight_layout()
    plt.savefig('ionic_conductivity.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Numerical examples
    print("=== Ionic Conductivity at SOFC Operating Temperature (800°C) ===")
    T_SOFC = 800 + 273.15
    for name, params in materials_ionic.items():
        sigma = ionic_conductivity(T_SOFC, params['sigma_0'], params['E_a'])
        print(f"{name:15s}: σ = {sigma:.4f} S/cm, E_a = {params['E_a']:.2f} eV")

# Execute
plot_ionic_conductivity()

Python Implementation: Temperature Dependence of YSZ Ionic Conductivity

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 6: Detailed Analysis of YSZ Conductivity
# ===================================

import numpy as np
import matplotlib.pyplot as plt

def ysz_conductivity_vs_doping():
    """
    Relationship between Y₂O₃ doping level and conductivity in YSZ
    """
    # Y₂O₃ doping level (mol%)
    doping_levels = np.array([3, 5, 8, 10, 12, 15])

    # Conductivity at 800°C (based on experimental data)
    sigma_800C = np.array([0.02, 0.06, 0.10, 0.08, 0.05, 0.03])

    # Activation energy
    E_a_values = np.array([1.1, 1.05, 1.0, 1.05, 1.15, 1.2])

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Doping level vs conductivity
    ax1 = axes[0]
    ax1.plot(doping_levels, sigma_800C, 'o-', linewidth=2, markersize=10, color='navy')
    ax1.axvline(x=8, color='red', linestyle='--', linewidth=2, label='Optimal: 8 mol%')

    ax1.set_xlabel('Y₂O₃ Content (mol%)', fontsize=12)
    ax1.set_ylabel('Ionic Conductivity σ (S/cm) @ 800°C', fontsize=12)
    ax1.set_title('YSZ Conductivity vs Doping Level', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # Right plot: Temperature dependence (for each doping level)
    ax2 = axes[1]
    T_range = np.linspace(600, 1200, 100) + 273.15

    for i, doping in enumerate([5, 8, 10]):
        idx = np.where(doping_levels == doping)[0][0]
        sigma_0_calc = sigma_800C[idx] / np.exp(-E_a_values[idx] / (8.617e-5 * (800 + 273.15)))
        sigma_T = ionic_conductivity(T_range, sigma_0_calc, E_a_values[idx])

        ax2.semilogy(T_range - 273.15, sigma_T, linewidth=2, label=f'{doping} mol% Y₂O₃')

    ax2.set_xlabel('Temperature (°C)', fontsize=12)
    ax2.set_ylabel('Ionic Conductivity σ (S/cm)', fontsize=12)
    ax2.set_title('YSZ Conductivity vs Temperature', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3, which='both')

    plt.tight_layout()
    plt.savefig('ysz_conductivity_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Numerical examples
    print("=== YSZ Doping Level Effects ===")
    for i, doping in enumerate(doping_levels):
        print(f"{doping:2d} mol% Y₂O₃: σ(800°C) = {sigma_800C[i]:.3f} S/cm, E_a = {E_a_values[i]:.2f} eV")

# Execute
ysz_conductivity_vs_doping()

4.5 Luminescent Material Ceramics

4.5.1 Luminescence Mechanism of Phosphors

Phosphors absorb excitation energy (electron beam, UV, blue LED) and emit visible light. Rare-earth ions (Ce³⁺, Eu²⁺, Eu³⁺, etc.) serving as luminescence centers are doped into host crystals (YAG, CaS, SrSiO₅, etc.).

Luminous efficiency (quantum efficiency) is:

\[ \eta = \frac{\text{Number of emitted photons}}{\text{Number of absorbed photons}} \times 100\% \]

4.5.2 YAG:Ce Phosphor for White LEDs

YAG:Ce (Y₃Al₅O₁₂:Ce³⁺) is excited by blue LEDs (450-470 nm) and emits yellow light (500-650 nm). White light is achieved through color mixing with the blue light.

Python Implementation: CIE Chromaticity Coordinate Calculation

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 7: CIE Chromaticity Coordinate Calculation
# ===================================

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simpson

def gaussian_emission_spectrum(wavelength, peak, fwhm, intensity=1.0):
    """
    Generate Gaussian emission spectrum

    Parameters:
    -----------
    wavelength : array
        Wavelength [nm]
    peak : float
        Peak wavelength [nm]
    fwhm : float
        Full width at half maximum (FWHM) [nm]
    intensity : float
        Intensity (arbitrary units)

    Returns:
    --------
    spectrum : array
        Emission intensity
    """
    sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))
    spectrum = intensity * np.exp(-((wavelength - peak)**2) / (2 * sigma**2))
    return spectrum


def cie_xyz_from_spectrum(wavelength, spectrum):
    """
    Calculate CIE XYZ values from emission spectrum (simplified version)

    Parameters:
    -----------
    wavelength : array
        Wavelength [nm]
    spectrum : array
        Emission intensity

    Returns:
    --------
    X, Y, Z : float
        CIE tristimulus values
    """
    # Simplified model of CIE 1931 standard observer (actual implementation uses standard data)
    # Gaussian approximation for simplicity

    # x̄(λ): Red cone sensitivity (peak 600 nm)
    x_bar = gaussian_emission_spectrum(wavelength, 600, 50, 1.0)

    # ȳ(λ): Green cone sensitivity (peak 550 nm)
    y_bar = gaussian_emission_spectrum(wavelength, 550, 50, 1.0)

    # z̄(λ): Blue cone sensitivity (peak 450 nm)
    z_bar = gaussian_emission_spectrum(wavelength, 450, 50, 1.0)

    # Integration
    X = simpson(spectrum * x_bar, wavelength)
    Y = simpson(spectrum * y_bar, wavelength)
    Z = simpson(spectrum * z_bar, wavelength)

    return X, Y, Z


def plot_yag_ce_emission():
    """
    YAG:Ce emission spectrum and CIE chromaticity coordinates
    """
    # Wavelength range
    wavelength = np.linspace(380, 780, 400)

    # Blue LED (excitation light)
    blue_led = gaussian_emission_spectrum(wavelength, 460, 20, 1.0)

    # YAG:Ce emission (yellow)
    yag_ce_emission = gaussian_emission_spectrum(wavelength, 560, 100, 0.8)

    # Combined spectrum (white light)
    white_light = blue_led + yag_ce_emission

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Emission spectrum
    ax1 = axes[0]
    ax1.fill_between(wavelength, 0, blue_led, alpha=0.5, color='blue', label='Blue LED (460 nm)')
    ax1.fill_between(wavelength, 0, yag_ce_emission, alpha=0.5, color='yellow', label='YAG:Ce (560 nm)')
    ax1.plot(wavelength, white_light, 'k-', linewidth=2, label='White Light (combined)')

    ax1.set_xlabel('Wavelength (nm)', fontsize=12)
    ax1.set_ylabel('Intensity (a.u.)', fontsize=12)
    ax1.set_title('YAG:Ce White LED Emission Spectrum', fontsize=14, fontweight='bold')
    ax1.set_xlim(380, 780)
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # Right plot: CIE 1931 chromaticity diagram (simplified version)
    ax2 = axes[1]

    # Calculate CIE chromaticity coordinates
    X_blue, Y_blue, Z_blue = cie_xyz_from_spectrum(wavelength, blue_led)
    x_blue = X_blue / (X_blue + Y_blue + Z_blue)
    y_blue = Y_blue / (X_blue + Y_blue + Z_blue)

    X_yellow, Y_yellow, Z_yellow = cie_xyz_from_spectrum(wavelength, yag_ce_emission)
    x_yellow = X_yellow / (X_yellow + Y_yellow + Z_yellow)
    y_yellow = Y_yellow / (X_yellow + Y_yellow + Z_yellow)

    X_white, Y_white, Z_white = cie_xyz_from_spectrum(wavelength, white_light)
    x_white = X_white / (X_white + Y_white + Z_white)
    y_white = Y_white / (X_white + Y_white + Z_white)

    # Plot
    ax2.scatter(x_blue, y_blue, s=200, c='blue', marker='o', edgecolors='black', linewidth=2, label='Blue LED')
    ax2.scatter(x_yellow, y_yellow, s=200, c='yellow', marker='s', edgecolors='black', linewidth=2, label='YAG:Ce')
    ax2.scatter(x_white, y_white, s=250, c='white', marker='*', edgecolors='black', linewidth=2, label='White Light', zorder=5)

    # Color temperature line (simplified)
    blackbody_x = np.array([0.45, 0.40, 0.35, 0.31])
    blackbody_y = np.array([0.41, 0.40, 0.38, 0.33])
    ax2.plot(blackbody_x, blackbody_y, 'k--', linewidth=1, alpha=0.5, label='Blackbody Locus')

    ax2.set_xlabel('CIE x', fontsize=12)
    ax2.set_ylabel('CIE y', fontsize=12)
    ax2.set_title('CIE 1931 Chromaticity Diagram', fontsize=14, fontweight='bold')
    ax2.set_xlim(0, 0.8)
    ax2.set_ylim(0, 0.9)
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('yag_ce_emission.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Numerical results
    print("=== CIE Chromaticity Coordinates ===")
    print(f"Blue LED:  (x, y) = ({x_blue:.3f}, {y_blue:.3f})")
    print(f"YAG:Ce:    (x, y) = ({x_yellow:.3f}, {y_yellow:.3f})")
    print(f"White Light: (x, y) = ({x_white:.3f}, {y_white:.3f})")
    print(f"\nColor Temperature: ~5000-6000 K (cool white)")

# Execute
plot_yag_ce_emission()

4.6 Magnetic Ceramics

4.6.1 Magnetism of Ferrites

Ferrites are oxide magnetic materials with the general formula MFe₂O₄ (M = Mn, Ni, Zn, Co). They have spinel structure and exhibit ferrimagnetism. Due to high electrical resistance (>10⁶ Ω·cm), they are suitable for high-frequency applications.

4.6.2 Key Property Parameters

• Saturation magnetization M_s: Maximum magnetization (T = Tesla)
• Coercivity H_c: Magnetic field to reverse magnetization (A/m)
• Permeability μ_r: Magnetic permeability (dimensionless)
• Curie temperature T_c: Ferromagnetic → paramagnetic transition temperature (°C)

Python Implementation: Ferrite Magnetization Property Simulation

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

# ===================================
# Example 8: Ferrite Magnetization Curves
# ===================================

import numpy as np
import matplotlib.pyplot as plt

def langevin_function(x):
    """
    Langevin function (paramagnetic magnetization)

    Parameters:
    -----------
    x : float or array
        Argument (μH/kT)

    Returns:
    --------
    L(x) : float or array
        Langevin function value
    """
    # Avoid divergence for small x
    with np.errstate(divide='ignore', invalid='ignore'):
        L = np.where(np.abs(x) < 1e-3, x / 3, 1 / np.tanh(x) - 1 / x)
    return L


def ferrite_magnetization_curve(H, M_s=0.35, H_c=15, a=1000):
    """
    Ferrite magnetization curve (simplified model)

    Parameters:
    -----------
    H : array
        Magnetic field [A/m]
    M_s : float
        Saturation magnetization [T]
    H_c : float
        Coercivity [A/m]
    a : float
        Shape parameter

    Returns:
    --------
    M : array
        Magnetization [T]
    """
    # tanh-based simplified model
    M = M_s * np.tanh((H - H_c) / a) if isinstance(H, (int, float)) else \
        M_s * np.tanh((H - H_c) / a)
    return M


def plot_ferrite_properties():
    """
    Visualize ferrite magnetic properties
    """
    # Magnetic field range
    H_range = np.linspace(-5000, 5000, 1000)

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left plot: Soft magnetic ferrite (MnZn)
    ax1 = axes[0]

    M_s_soft = 0.35
    H_c_soft = 15
    M_soft = M_s_soft * np.tanh(H_range / 500)

    ax1.plot(H_range, M_soft, linewidth=2, color='navy')
    ax1.axhline(y=M_s_soft, color='red', linestyle='--', linewidth=1, alpha=0.5, label=f'M_s = {M_s_soft} T')
    ax1.axhline(y=-M_s_soft, color='red', linestyle='--', linewidth=1, alpha=0.5)
    ax1.axvline(x=0, color='black', linestyle='-', linewidth=0.5)
    ax1.axhline(y=0, color='black', linestyle='-', linewidth=0.5)

    ax1.set_xlabel('Magnetic Field H (A/m)', fontsize=12)
    ax1.set_ylabel('Magnetization M (T)', fontsize=12)
    ax1.set_title('Soft Ferrite (MnZn) M-H Curve', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # Right plot: Hard magnetic ferrite (BaFe₁₂O₁₉)
    ax2 = axes[1]

    M_s_hard = 0.48
    H_c_hard = 250000  # 250 kA/m
    H_hard = np.linspace(-400000, 400000, 1000)
    M_hard = M_s_hard * np.tanh((H_hard - H_c_hard) / 50000)

    # Hysteresis loop approximation
    M_hard_up = M_s_hard * np.tanh((H_hard + H_c_hard) / 50000)
    M_hard_down = M_s_hard * np.tanh((H_hard - H_c_hard) / 50000)

    ax2.plot(H_hard / 1000, M_hard_up, linewidth=2, color='blue', label='Ascending')
    ax2.plot(H_hard / 1000, M_hard_down, linewidth=2, color='red', label='Descending')
    ax2.axvline(x=H_c_hard / 1000, color='green', linestyle='--', linewidth=2, label=f'H_c = {H_c_hard/1000:.0f} kA/m')
    ax2.axvline(x=-H_c_hard / 1000, color='green', linestyle='--', linewidth=2)

    ax2.set_xlabel('Magnetic Field H (kA/m)', fontsize=12)
    ax2.set_ylabel('Magnetization M (T)', fontsize=12)
    ax2.set_title('Hard Ferrite (BaFe₁₂O₁₉) Hysteresis Loop', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('ferrite_magnetization.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Numerical examples
    print("=== Ferrite Magnetic Properties ===")
    print("\nSoft Ferrite (MnZn):")
    print(f"  M_s = {M_s_soft} T, H_c = {H_c_soft} A/m")
    print(f"  Application: High permeability (μ_r ~ 2000-5000), transformer cores")

    print("\nHard Ferrite (BaFe₁₂O₁₉):")
    print(f"  M_s = {M_s_hard} T, H_c = {H_c_hard/1000} kA/m")
    print(f"  Application: Permanent magnets, motors, speakers")

# Execute
plot_ferrite_properties()

Exercises

alpha = 2.5e-40
N = 5e28
epsilon_r = clausius_mossotti_permittivity(alpha, N)
print(f"εᵣ = {epsilon_r:.1f}")
# Output: εᵣ ≈ 80

C_target = 1e-6  # 1 μF
epsilon_r = 3000
A = (1e-3)**2
d = 5e-6
epsilon_0 = 8.854e-12

n = (C_target * d) / (epsilon_0 * epsilon_r * A)
print(f"Required number of layers n = {n:.0f}")
# Output: n ≈ 377 layers

d_33 = 593e-12
E = 1e6
thickness = 1e-3

strain = piezoelectric_strain(d_33, E)
displacement = strain * thickness
print(f"Displacement = {displacement*1e9:.1f} nm")
# Output: Displacement ≈ 593 nm

T1 = 700 + 273.15
T2 = 900 + 273.15
sigma1 = 0.05
sigma2 = 0.15
k_B = 8.617e-5

# ln(σ₂/σ₁) = -E_a/k_B × (1/T₂ - 1/T₁)
E_a = -k_B * np.log(sigma2 / sigma1) / (1/T2 - 1/T1)
print(f"E_a = {E_a:.2f} eV")
# Output: E_a ≈ 0.95 eV

# To increase color temperature (cool white), increase blue component
# Reduce YAG:Ce emission intensity or make YAG:Ce layer thinner

# 4000 K: Blue/Yellow = 1.0
# 6500 K: Blue/Yellow = 1.5

print("To increase color temperature:")
print("1. Make YAG:Ce phosphor layer thinner (reduce yellow)")
print("2. Reduce YAG:Ce concentration")
print("3. Increase blue LED output")

x_Zr = np.linspace(0, 1, 100)
x_MPB = 0.52
d_33 = 200 + 400 * np.exp(-((x_Zr - x_MPB) / 0.1)**2)

max_idx = np.argmax(d_33)
print(f"Optimal composition: x_Zr = {x_Zr[max_idx]:.2f}")
print(f"Maximum d₃₃ = {d_33[max_idx]:.0f} pC/N")
# Output: x_Zr = 0.52, d₃₃ ≈ 600 pC/N

print("Causes of frequency dependence:")
print("1. Eddy current loss (low resistivity)")
print("2. Delayed domain wall motion (domain walls cannot follow at high frequency)")
print("3. Spin rotation resonance (ferrimagnetic resonance frequency)")
print("\nCountermeasure: Switch to NiZn ferrite (high resistance)")

sigma = 0.1  # S/cm = 10 S/m
rho = 1 / sigma  # Ω·m
A = 10e-4  # m²
R_target = 0.1  # Ω

L_max = R_target * A / rho
print(f"Maximum thickness L = {L_max*1e6:.0f} μm")
# Output: L ≈ 100 μm
# In actual SOFCs, 10-50 μm is typical

d_33 = 300e-12
F = 100  # N
A = 10e-4  # m²
f = 10  # Hz
epsilon_r = 1300
epsilon_0 = 8.854e-12

# Generated charge Q = d₃₃ × F
Q = d_33 * F
print(f"Generated charge Q = {Q*1e9:.1f} nC")

# Capacitance C = ε₀εᵣA/d (assuming 1 mm thickness)
d = 1e-3
C = epsilon_0 * epsilon_r * A / d
V = Q / C
print(f"Voltage V = {V:.1f} V")

# Power P = 0.5 × C × V² × f
P = 0.5 * C * V**2 * f
print(f"Power P = {P*1e6:.1f} μW")
# Example output: ~10-100 μW (practical level requires mW)

print("Design strategy:")
print("1. YAG:Ce layer (yellow): 550 nm emission with blue LED excitation")
print("2. CaAlSiN₃:Eu layer (red): 630 nm emission with blue excitation")
print("3. Thickness ratio: YAG:Ce = 80%, red phosphor = 20%")
print("4. Emission intensity ratio: Blue:Yellow:Red = 30:50:20")
print("5. Adjust CIE chromaticity coordinates to (0.33, 0.33)")
print("6. Confirm color rendering index Ra > 90 (R9 red component is critical)")
print("\nOptimization algorithm:")
print("- Monte Carlo optimization with thickness and concentration as parameters")
print("- Objective function: Maximize Ra & color temperature constraint 5000-6500K")

References

1. Moulson, A.J., Herbert, J.M. (2003). Electroceramics: Materials, Properties, Applications. Wiley, pp. 1-85, 155-210, 340-395.
2. Jaffe, B., Cook, W.R., Jaffe, H. (1971). Piezoelectric Ceramics. Academic Press, pp. 50-135.
3. Tuller, H.L. (2000). Ionic conduction in nanocrystalline materials. Solid State Ionics, 131, 15-68.
4. Ronda, C.R. (2008). Luminescence: From Theory to Applications. Wiley-VCH, pp. 120-180.
5. Blasse, G., Grabmaier, B.C. (1994). Luminescent Materials. Springer, pp. 1-65.
6. Goldman, A. (2006). Modern Ferrite Technology. Springer, pp. 30-95.
7. Setter, N. (2002). Piezoelectric Materials in Devices. EPFL Swiss Federal Institute of Technology, pp. 1-50, 200-250.