The Mullins–Herring (MH) Model is a theoretical framework used in surface science and materials physics to describe surface diffusion and the smoothing of interfaces. Developed to analyze the kinetics of surface morphological changes, the model applies to scenarios such as crystal growth, thin film deposition, and annealing processes. It operates by considering how atoms migrate across a surface to reduce surface free energy, predicting the evolution of features like ripples and mounds over time. The MH Model emphasizes the importance of curvature-driven diffusion, explaining how atoms move from high-curvature regions to lower-curvature areas, leading to the gradual flattening of surfaces.

🌵MATLAB snippet

The Mullins–Herring (MH) model is a higher-order SPDE often used to describe surface diffusion during growth processes. It models the surface evolution with fourth-order spatial derivatives, which is a step up in complexity from the KPZ equation. The MH equation can be expressed as:

$$ \frac{\partial h(x, t)}{\partial t} = -\nu \nabla^4 h(x, t) + \eta(x, t), $$

where:

• ( $h(x, t)$ ) is the height function,
• ( $\nu$ ) is the diffusion coefficient,
• ( $\nabla^4$ ) represents the biharmonic operator (fourth-order spatial derivative),
• ( $\eta(x, t)$ ) is the Gaussian white noise.

MATLAB Implementation Steps:

1. Discretize the spatial domain: Create a grid to represent ( $h(x, t)$ ).
2. Finite difference scheme: Use a fourth-order central difference to approximate ( $\nabla^4 h(x, t)$ ).
3. Time integration: Implement a simple explicit scheme like Euler-Maruyama.

Here's a basic MATLAB script to simulate the Mullins–Herring model:

 % Parameters
 N = 100;               % Number of spatial points
 L = 10;                % Length of the domain
 dx = L / N;            % Spatial step size
 dt = 0.01;             % Time step size
 nu = 1.0;              % Diffusion coefficient
 total_time = 1.0;      % Total simulation time
 noise_strength = 0.1;  % Strength of the noise
 ​
 % Grid setup
 x = linspace(0, L, N);
 h = zeros(1, N);       % Initial height function
 ​
 % Noise generator
 rng(0); % Set random seed for reproducibility
 ​
 % Pre-compute indices for finite differences
 indices = 1:N;
 ip1 = mod(indices, N) + 1; % i+1 with periodic boundary
 im1 = mod(indices - 2, N) + 1; % i-1 with periodic boundary
 ip2 = mod(indices + 1, N) + 1; % i+2 with periodic boundary
 im2 = mod(indices - 3, N) + 1; % i-2 with periodic boundary
 ​
 % Time evolution loop
 num_steps = total_time / dt;
 for t = 1:num_steps
     % Compute the biharmonic term using fourth-order finite differences
     laplacian2 = (-h(im2) + 4*h(im1) - 6*h(indices) + 4*h(ip1) - h(ip2)) / (dx^4);

     % Add noise term
     noise = noise_strength * sqrt(dt) * randn(1, N);

     % Update rule (Euler-Maruyama step)
     h = h + dt * (-nu * laplacian2) + noise;

     % Optional: Plot every few steps to visualize
     if mod(t, 100) == 0
         plot(x, h, 'b');
         xlabel('x');
         ylabel('Height h(x)');
         title(['Time step: ', num2str(t * dt)]);
         drawnow;
     end
 end

Explanation:

• Fourth-order finite difference: The biharmonic operator ( $\nabla^4$ ) is approximated with a combination of shifted values around each grid point.
• Noise: Gaussian random noise is added at each time step, scaled by ( $\sqrt{\text{dt}}$ ) to maintain the correct stochastic behavior.
• Periodic boundary conditions: The use of mod ensures that the grid wraps around at the boundaries.

Visualization: