This program was contributed by Daniel Garcia-Sanchez.
Note
As a prerequisite of this program, you need to have HDF5, complex PETSc, and the p4est libraries installed. The installation of deal.II together with these additional libraries is described in the README file.
Introduction
A phononic crystal is a periodic nanostructure that modifies the motion of mechanical vibrations or phonons. Phononic structures can be used to disperse, route and confine mechanical vibrations. These structures have potential applications in quantum information and have been used to study macroscopic quantum phenomena. Phononic crystals are usually fabricated in cleanrooms.
In this tutorial we show how to a design a phononic superlattice cavity which is a particular type of phononic crystal that can be used to confine mechanical vibrations. A phononic superlattice cavity is formed by two Distributed Bragg Reflector, mirrors and a \(\lambda/2\) cavity where \(\lambda\) is the acoustic wavelength. Acoustic DBRs are periodic structures where a set of bilayer stacks with contrasting physical properties (sound velocity index) is repeated \(N\) times. Superlattice cavities are usually grown on a Gallium Arsenide wafer by Molecular Beam Epitaxy. The bilayers correspond to GaAs/AlAs mirror pairs. As shown below, the thickness of the mirror layers (brown and green) is \(\lambda/4\) and the thickness of the cavity (blue) is \(\lambda/2\).
In this tutorial we calculate the band gap and the mechanical resonance of a phononic superlattice cavity but the code presented here can be easily used to design and calculate other types of phononic crystals.
The device is a waveguide in which the wave goes from left to right. The simulations of this tutorial are done in 2D, but the code is dimension independent and can be easily used with 3D simulations. The waveguide width is equal to the \(y\) dimension of the domain and the waveguide length is equal to the \(x\) dimension of the domain. There are two regimes that depend on the waveguide width:
Single mode: In this case the width of the structure is much smaller than the wavelength. This case can be solved either with FEM (the approach that we take here) or with a simple semi-analytical 1D transfer matrix formalism.
Multimode: In this case the width of the structure is larger than the wavelength. This case can be solved using FEM or with a scattering matrix formalism. Although we do not study this case in this tutorial, it is very easy to reach the multimode regime by increasing the parameter waveguide width (dimension_y in the jupyter notebook).
The simulations of this tutorial are performed in the frequency domain. To calculate the transmission spectrum, we use a procedure that is commonly used in time domain FDTD simulations. A pulse at a certain frequency is generated on the left side of the structure and the transmitted energy is measured on the right side of the structure. The simulation is run twice. First, we run the simulation with the phononic structure and measure the transmitted energy:
Then, we run the simulation without the phononic structure and measure the transmitted energy. We use the simulation without the structure for the calibration:
The transmission coefficient corresponds to the energy of the first simulation divided by the calibration energy. We repeat this procedure for each frequency step.
Elastic equations
What we want to simulate here is the transmission of elastic waves. Consequently, the right description of the problem uses the elastic equations, which in the time domain are given by
A perfectly matched layer (PML) can be used to truncate the solution at the boundaries. A PML is a transformation that results in a complex coordinate stretching.
Instead of a time domain approach, this tutorial program converts the equations above into the frequency domain by performing a Fourier transform with regard to the time variable. The elastic equations in the frequency domain then read as follows
where the coefficients \(s_i = 1+is_i'(x,y,z)\) account for the absorption. There are 3 \(s_i\) coefficients in 3D and 2 in 2D. The imaginary par of \(s_i\) is equal to zero outside the PML. The PMLs are reflectionless only for the exact wave equations. When the set of equations is discretized the PML is no longer reflectionless. The reflections can be made arbitrarily small as long as the medium is slowly varying, see the adiabatic theorem. In the code a quadratic turn-on of the PML has been used. A linear and cubic turn-on is also known to work. These equations can be expanded into
where summation over repeated indices (here \(n\), as well as \(k\) and \(l\)) is as always implied. Note that the strain is no longer symmetric after applying the complex coordinate stretching of the PML. This set of equations can be written as
The same as the strain, the stress tensor is not symmetric inside the PML ( \(s_j\neq 0\)). Indeed the fields inside the PML are not physical. It is useful to introduce the tensors \(\alpha_{mnkl}\) and \(\beta_{mnkl}\).
It is this set of equations we want to solve for a set of frequencies \(\omega\) in order to compute the transmission coefficient as function of frequency. The linear system becomes
In this tutorial we use a python jupyter notebook to set up the parameters and run the simulation. First we create a HDF5 file where we store the parameters and the results of the simulation.
Each of the simulations (displacement and calibration) is stored in a separate HDF5 group:
import numpy as np
import h5py
import matplotlib.pyplot as plt
import subprocess
import scipy.constants as constants
import scipy.optimize
# This considerably reduces the size of the svg data
plt.rcParams['svg.fonttype'] = 'none'
h5_file = h5py.File('results.h5', 'w')
data = h5_file.create_group('data')
displacement = data.create_group('displacement')
calibration = data.create_group('calibration')
# Set the parameters
for group in [displacement, calibration]:
# Dimensions of the domain
# The waveguide length is equal to dimension_x
group.attrs['dimension_x'] = 2e-5
# The waveguide width is equal to dimension_y
group.attrs['dimension_y'] = 2e-8
# Position of the probe that we use to measure the flux
The variable data is the HDF5::Group in which all the simulation results will be stored. Note that the variables RightHandSide::data, PML::data, Rho::data and Parameters::data point to the same group of the HDF5 file. When a HDF5::Group is copied, it will point to the same group of the HDF5 file.
The simulation parameters are stored in data as HDF5 attributes. The following attributes are defined in the jupyter notebook, stored in data as HDF5 attributes and then read by the constructor.
The calculation of the mass and stiffness matrices is very expensive. These matrices are the same for all the frequency steps. The right hand side vector is also the same for all the frequency steps. We use this class to store these objects and re-use them at each frequency step. Note that here we don't store the assembled mass and stiffness matrices and right hand sides, but instead the data for a single cell. QuadratureCache class is very similar to the PointHistory class that has been used in step-18.
template <int dim>
class QuadratureCache
{
public:
QuadratureCache(constunsignedint dofs_per_cell);
private:
unsignedint dofs_per_cell;
public:
We store the mass and stiffness matrices in the variables mass_coefficient and stiffness_coefficient. We store as well the right_hand_side and JxW values which are going to be the same for all the frequency steps.
This function returns the stiffness tensor of the material. For the sake of simplicity we consider the stiffness to be isotropic and homogeneous; only the density \(\rho\) depends on the position. As we have previously shown in step-8, if the stiffness is isotropic and homogeneous, the stiffness coefficients \(c_{ijkl}\) can be expressed as a function of the two coefficients \(\lambda\) and \(\mu\). The coefficient tensor reduces to
where \(a\) is the maximum amplitude that takes the force and \(\sigma_x\) and \(\sigma_y\) are the standard deviations for the \(x\) and \(y\) components. Note that the pulse has been cropped to \(x_\textrm{min}<x<x_\textrm{max}\) and \(y_\textrm{min} <y<y_\textrm{max}\).
As before, the constructor reads all the parameters from the HDF5::Groupdata using the HDF5::Group::get_attribute() function. As we have discussed, a quadratic turn-on of the PML has been defined in the jupyter notebook. It is possible to use a linear, cubic or another power degree by changing the parameter pml_coeff_degree. The parameters pml_x and pml_y can be used to turn on and off the x and y PMLs.
This class is used to define the mass density. As we have explaine before, a phononic superlattice cavity is formed by two Distributed Reflector, mirrors and a \(\lambda/2\) cavity where \(\lambda\) is the acoustic wavelength. Acoustic DBRs are periodic structures where a set of bilayer stacks with contrasting physical properties (sound velocity index) is repeated \(N\) times. The change of in the wave velocity is generated by alternating layers with different density.
where \(K_e\) is the effective elastic constant and \(\rho\) the density. Here we consider the case in which the waveguide width is much smaller than the wavelength. In this case it can be shown that for the two dimensional case
\[ K_e = 4\mu\frac{\lambda +\mu}{\lambda+2\mu} \]
and for the three dimensional case \(K_e\) is equal to the Young's modulus.
This is very similar to the constructor of step-40. In addition we create the HDF5 datasets frequency_dataset, position_dataset and displacement. Note the use of the template keyword for the creation of the HDF5 datasets. It is a C++ requirement to use the template keyword in order to treat create_dataset as a dependent template name.
There is nothing new in this function, the only difference with step-40 is that we don't have to apply boundary conditions because we use the PMLs to truncate the domain.
This function is also very similar to step-40, though there are notable differences. We assemble the system for each frequency/omega step. In the first step we set calculate_quadrature_data = True and we calculate the mass and stiffness matrices and the right hand side vector. In the subsequent steps we will use that data to accelerate the calculation.
We calculate the stiffness tensor for the \(\lambda\) and \(\mu\) that have been defined in the jupyter notebook. Note that contrary to \(\rho\) the stiffness is constant among for the whole domain.
We have to calculate the values of the right hand side, rho and the PML only if we are going to calculate the mass and the stiffness matrices. Otherwise we can skip this calculation which considerably reduces the total calculation time.
We have done this in step-18. Get a pointer to the quadrature cache data local to the present cell, and, as a defensive measure, make sure that this pointer is within the bounds of the global array:
Here we calculate the stiffness matrix. Note that the stiffness matrix is not symmetric because of the PMLs. We use the gradient function (see the documentation) which is a Tensor<2,dim>. The matrix \(G_{ij}\) consists of entries
Note the position of the indices \(i\) and \(j\) and the notation that we use in this tutorial: \(\partial_j\phi_i\). As the stiffness tensor is not symmetric, it is very easy to make a mistake.
We save the value of the stiffness matrix in quadrature_data
quadrature_data.stiffness_coefficient[i][j] =
stiffness_coefficient;
}
and the value of the right hand side in quadrature_data.
quadrature_data.right_hand_side[i] =
phi_i * force * fe_values.JxW(q);
}
}
We loop again over the degrees of freedom of the cells to calculate the system matrix. These loops are really quick because we have already calculated the stiffness and mass matrices, only the value of \(\omega\) changes.
This is even more simple than in step-40. We use the parallel direct solver MUMPS which requires less options than an iterative solver. The drawback is that it does not scale very well. It is not straightforward to solve the Helmholtz equation with an iterative solver. The shifted Laplacian multigrid method is a well known approach to precondition this system, but this is beyond the scope of this tutorial.
The vector coordinates contains the coordinates in the HDF5 file of the points of the probe that are located in locally owned cells. The vector displacement_data contains the value of the displacement at these points.
We write the displacement data in the HDF5 file. The call HDF5::DataSet::write_selection() is MPI collective which means that all the processes have to participate.
Therefore even if the process has no data to write it has to participate in the collective call. For this we can use HDF5::DataSet::write_none(). Note that we have to specify the data type, in this case std::complex<double>.
else
{
displacement.write_none<std::complex<double>>();
}
If the variable save_vtu_files in the input file equals True then all the data will be saved as vtu. The procedure to write vtu files has been described in step-40.
And on the cells that we are not interested in, set the respective value to a bogus value in order to make sure that if we were somehow wrong about our assumption we would find out by looking at the graphical output:
else
{
for (unsignedint dim_idx = 0; dim_idx < dim; ++dim_idx)
This function writes the datasets that have not already been written.
template <int dim>
void ElasticWave<dim>::output_results()
{
The vectors frequency and position are the same for all the processes. Therefore any of the processes can write the corresponding datasets. Because the call HDF5::DataSet::write is MPI collective, the rest of the processes will have to call HDF5::DataSet::write_none.
We use this function at the beginning of our computations to set up initial values of the cache variables. This function has been described in step-18. There are no differences with the function of step-18.
For clarity we divide the function run of step-40 into the functions run and frequency_sweep. In the function frequency_sweep we place the iteration over the frequency vector.
template <int dim>
void ElasticWave<dim>::frequency_sweep()
{
for (unsignedint frequency_idx = 0;
frequency_idx < parameters.nb_frequency_points;
++frequency_idx)
{
pcout << parameters.simulation_name + " frequency idx: "
In the first frequency step we calculate the mass and stiffness matrices and the right hand side. In the subsequent frequency steps we will use those values. This improves considerably the calculation time.
A phononic cavity is characterized by the resonance frequency and the the quality factor. The quality factor is equal to the ratio between the stored energy in the resonator and the energy dissipated energy per cycle, which is approximately equivalent to the ratio between the resonance frequency and the full width at half maximum (FWHM). The FWHM is equal to the bandwidth over which the power of vibration is greater than half the power at the resonant frequency.
where \(f_r = \frac{\omega_r}{2\pi}\) is the resonance frequency and \(\Gamma=\frac{\omega_r}{Q}\) is the dissipation rate. We used the previous equation in the jupyter notebook to fit the mechanical resonance.
Given the values we have chosen for the parameters, one could estimate the resonance frequency analytically. Indeed, this is then confirmed by what we get in this program: the phononic superlattice cavity exhibits a mechanical resonance at 20GHz and a quality factor of 5046. The following images show the transmission amplitude and phase as a function of frequency in the vicinity of the resonance frequency:
The images above suggest that the periodic structure has its intended effect: It really only lets waves of a very specific frequency pass through, whereas all other waves are reflected. This is of course precisely what one builds these sorts of devices for. But it is not quite this easy. In practice, there is really only a "band gap", i.e., the device blocks waves other than the desired one at 20GHz only within a certain frequency range. Indeed, to find out how large this "gap" is within which waves are blocked, we can extend the frequency range to 16 GHz through the appropriate parameters in the input file. We then obtain the following image:
What this image suggests is that in the range of around 18 to around 22 GHz, really only the waves with a frequency of 20 GHz are allowed to pass through, but beyond this range, there are plenty of other frequencies that can pass through the device.
Mode profile
We can inspect the mode profile with Paraview or VisIt. As we have discussed, at resonance all the mechanical energy is transmitted and the amplitude of motion is amplified inside the cavity. It can be observed that the PMLs are quite effective to truncate the solution. The following image shows the mode profile at resonance:
On the other hand, out of resonance all the mechanical energy is reflected. The following image shows the profile at 19.75 GHz. Note the interference between the force pulse and the reflected wave at the position \(x=-8\mu\textrm{m}\).
Experimental applications
Phononic superlattice cavities find application in quantum optomechanics. Here we have presented the simulation of a 2D superlattice cavity, but this code can be used as well to simulate "real world" 3D devices such as micropillar superlattice cavities, which are promising candidates to study macroscopic quantum phenomena. The 20GHz mode of a micropillar superlattice cavity is essentially a mechanical harmonic oscillator that is very well isolated from the environment. If the device is cooled down to 20mK in a dilution fridge, the mode would then become a macroscopic quantum harmonic oscillator.
1.8.17