7.2. Introduction to piezoelectric properties with lithium niobate

Continuing the study of bulk material elastic properties from the previous tutorial, let’s take a look at the orientational dependence of lithium niobate, an increasingly important material in SBS.

The numerical values for the \(\mathrm{LiNbO}_3\) material properties are taken from Rodrigues et al, JOSA B 40, D56 (2023).

%load_ext autoreload
%autoreload 3

import sys
#import os
import time
import numpy as np
from IPython.display import Image, display, HTML

sys.path.append("../backend")
import materials

def img_single(fn, cap='', width=300):
    cache_buster=int(time.time() * 1000)
    #<div style="display: flex; justify-content: center; align-items: flex-start;">
    #</div>

    html_string = f"""
        <div style="margin-right: 10px;">
            <img src="{fn}?{cache_buster}" width="{width}px">
            <p style="text-align: center;">{cap}</p>
        </div>"""

    display(HTML(html_string))

def img_pair(fn1, fn2, cap1='', cap2='', width1=300, width2=300):
    cache_buster=int(time.time() * 1000)
    html_string = f"""
    <div style="display: flex; justify-content: center; align-items: flex-start;">
        <div style="margin-right: 10px;">
            <img src="{fn1}?{cache_buster}" width="{width1}px">
            <p style="text-align: center;">{cap1}</p>
        </div>
        <div>
            <img src="{fn2}?{cache_buster}" width="{width2}px">
            <p style="text-align: center;">{cap2}</p>
        </div>
    </div>"""

    display(HTML(html_string))

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

7.2.1. Stiffness tensor under crystal rotations

7.2.1.1. The default orientation: \(y\)-cut

To review, in NumBAT, the laboratory axes \(x\), \(y\) and \(z\) are fixed. Propagation always occurs along the \(z\) direction and the vertical direction out of the substrate is always along \(y\).

The crystal axes of a material, denoted \(X\), \(Y\) and \(Z\), or \(\hat{c}_x\), \(\hat{c}_y\), \(\hat{c}_z\) can be rotated.

The material is defined in the definition json files with the crystal and laboratory axes aligned, and is thus in \(y\)-cut form: the crystal \(\hat{c}_y\) axis points up out of the waveguide, with the crystal symmetry axis \(\hat{c}_z\) lying along the propagation direction \(z\).

We’ll look at this for a material file describing lithium niobate. (Remember that NumBAT, the same material can have multiple definition files to allow exact comparison with literature references.)

mat_LiNbO3 = materials.make_material("LiNbO3_2023_Rodrigues")
mat_LiNbO3_y = mat_LiNbO3  # for naming convenience later

Here is the default orientation as just described with the crystal and laboratory axes aligned.

pref='tmp_LiNbO3'
fnimg = mat_LiNbO3.make_crystal_axes_plot(pref+'-ycut')
img_single(fnimg, width=300)

The crystal properties show that for this orientation the two shear modes for propagation along \(z\) are degenerate:

print(mat_LiNbO3.elastic_properties())

Elastic properties of material LiNbO3_2023_Rodrigues
  Density:        4650.000 kg/m^3
  Ref. index:     2.2100+0.0000j
  Crystal class:  Trigonal
  Crystal group:  no sym
  Stiffness c_IJ:
 Stiffness c_IJ, unit: GPa.
   Voigt 4-tensor:
    [[198.8300  54.6400  68.2300   7.8300   0.0000   0.0000]
     [ 54.6400 198.8300  68.2300  -7.8300   0.0000   0.0000]
     [ 68.2300  68.2300 235.7100   0.0000   0.0000   0.0000]
     [  7.8300  -7.8300   0.0000  59.8600   0.0000   0.0000]
     [  0.0000   0.0000   0.0000   0.0000  59.8600   7.8300]
     [  0.0000   0.0000   0.0000   0.0000   7.8300  72.0950]]

  Wave mode 1: v_p=7.1197 km/s,  |v_g|=7.1197 km/s,  u_j=[ 0.0000  0.0000  1.0000],  v_g=[ 0.0000  0.0000  7.1197] km/s
  Wave mode 2: v_p=3.5879 km/s,  |v_g|=3.6185 km/s,  u_j=[ 1.0000  0.0000  0.0000],  v_g=[ 0.0000  0.4693  3.5879] km/s
  Wave mode 3: v_p=3.5879 km/s,  |v_g|=3.6185 km/s,  u_j=[ 0.0000  1.0000  0.0000],  v_g=[ 0.0000 -0.4693  3.5879] km/s

This is not true if the wave vector direction in the \(x\)-\(z\) plane is varied.

By solving the Christoffel equation as the wave vector rotates in the \(x\)-\(z\) plane we can map out curves of the phase velocity, inverse phase velocity (or “slowness”), and the magnitude of the group velocity \(|v_g(\kappa)|\).

The faint circular lines mark radial velocities at 1 km/s intervals.

imgdisp_y = mat_LiNbO3.plot_bulk_dispersion(pref+'-ycut')
img_single(imgdisp_y, width=400)

For this material, it is quite hard to see the difference between the phase and group velocity plots. That is more obvious in some other materials.

These plots can also be generated without the polarisation state markers:

imgdisp_yb = mat_LiNbO3.plot_bulk_dispersion(pref, show_poln=False)
img_single(imgdisp_yb, width=400)

7.2.1.2. The \(z\)-cut orientation

Now let’s rotate the starting crystal so the primary crystal symmetry axis \(\hat{c}_z\) points upwards along \(\hat{y}\), with the \(\hat{c}_y\) axis along \(-\hat{z}\). This corresponds to Fig. 2a in Rodrigues et al.

mat_LiNbO3_z = mat_LiNbO3.copy()
mat_LiNbO3_z.set_orientation('z-cut')

fimg = mat_LiNbO3_z.make_crystal_axes_plot(pref+'-zcut')
img_single(fimg, width=300)

With the 6-fold symmetry axis \(Z\) pointing up along \(y\), the dispersion cuts in the \(x\)-\(z\) plane now display the full 6-fold symmetry associated with the trigonal structure of lithium niobate:

imgdisp_z = mat_LiNbO3_z.plot_bulk_dispersion(pref+'-zcut')
img_single(imgdisp_z,width=500)

7.2.1.3. The \(x\)-cut orientation

For completeness, we can look at the \(x\)-cut properties:

mat_LiNbO3_x =  mat_LiNbO3.copy()
mat_LiNbO3_x.set_orientation('x-cut')

print(mat_LiNbO3_x.elastic_properties())

Elastic properties of material LiNbO3_2023_Rodrigues
  Density:        4650.000 kg/m^3
  Ref. index:     2.2100+0.0000j
  Crystal class:  Trigonal
  Crystal group:  no sym
  Stiffness c_IJ:
 Stiffness c_IJ, unit: GPa.
   Voigt 4-tensor:
    [[198.8300  54.6400  68.2300  -0.0000   7.8300   0.0000]
     [ 54.6400 198.8300  68.2300   0.0000  -7.8300  -0.0000]
     [ 68.2300  68.2300 235.7100   0.0000   0.0000   0.0000]
     [ -0.0000   0.0000   0.0000  59.8600   0.0000  -7.8300]
     [  7.8300  -7.8300   0.0000   0.0000  59.8600  -0.0000]
     [  0.0000  -0.0000   0.0000  -7.8300  -0.0000  72.0950]]

  Wave mode 1: v_p=7.1197 km/s,  |v_g|=7.1197 km/s,  u_j=[ 0.0000  0.0000  1.0000],  v_g=[ 0.0000  0.0000  7.1197] km/s
  Wave mode 2: v_p=3.5879 km/s,  |v_g|=3.6185 km/s,  u_j=[ 1.0000  0.0000  0.0000],  v_g=[ 0.4693 -0.0000  3.5879] km/s
  Wave mode 3: v_p=3.5879 km/s,  |v_g|=3.6185 km/s,  u_j=[ 0.0000  1.0000  0.0000],  v_g=[-0.4693  0.0000  3.5879] km/s

fnimg = mat_LiNbO3_x.make_crystal_axes_plot(pref+'-xcut')
img_single(fnimg,width=300)

The distinction between the \(x\)-cut and \(y\)-cut dispersion is subtle, and it helpful to plot them side by side.

imgdisp_x = mat_LiNbO3_x.plot_bulk_dispersion(pref+'-xcut')
img_pair(imgdisp_y, imgdisp_x, 'Y cut', 'X cut', 400, 400)

Y cut

X cut

7.2.2. Photoelastic response

Now let’s consider the influence of crystal orientation on the photoelastic coupling.

The photoelastic tensor elements \(p_{ijkl}\) also change in value as the crystal is rotated. Once again, we can plot the coefficients in the lab frame coordinates, as a crystal of given orientation is imagined to be rotated around the \(y\) axis.

The desired elements of the tensor are specified in the Voigt notation \(p_{IJ}\).

As with the results above, the \(Z\)-cut crystal exhibits the 3-fold trigonal symmetry, while the other two cuts have more complex behaviour.

mat_LiNbO3_z.plot_photoelastic_IJ(pref, ("11","12", "13", "14", "15"))

mat_LiNbO3_x.plot_photoelastic_IJ(pref, ("11","12", "13", "14", "15"))

mat_LiNbO3_y.plot_photoelastic_IJ(pref, ("11","12", "13", "14", "15"))

7.3. Piezolectric properties of lithium niobate

Now let’s turn to the piezoelectric response of lithium niobate.

This version of lithium niobate has additional piezoelectric quantities defined: the strain piezo coefficient \(d_{iJ}\) and the stress piezo coefficient \(e_{iJ}\). These are \(3\times 6\) tensors with mixed ordinary and Voigt indices.

7.3.1. Piezo slowness curves for x-cut

There are a number of piezoelectric phenomena that involve coupling of the optical field and a material-related polarisation field. Perhaps the most basic effect is a “piezo-stiffening” of the stiffness tensor (a genuinely confusing terminology). This effect is discussed in Auld volume 1, section 8.F.

To good approximation this is an adjustment to the optical properties only in which the “raw” stiffness tensor \(c_{IJ}\) is modified by the piezo influence.

It is described by a change to the Christoffel equation so that the bulk optical modes are altered. The revised Christoffel equation has the form

\[q^2 \Gamma_{ij} v_j = \rho \Omega^2 v_i,\]

where

\[\Gamma_{ij} = l_{iK} c_{KL} l_{Lj},\]

and the matrix \(l_{Lj}\) indexed by Voigt (\(L\)) and Cartesian indices (\(j\)) is composed of the wavevector components $:nbsphinx-math:vec `:nbsphinx-math:kappa `$:

\[\begin{split}l_{Lj} = \begin{bmatrix} \kappa_x & 0 & 0 \\ 0 & \kappa_y & 0 \\ 0 & 0& \kappa_z \\ 0 & \kappa_z & \kappa_y \\ \kappa_z & 0 & \kappa_x \\ \kappa_y & \kappa_x & 0 \end{bmatrix}.\end{split}\]

The piezoelastically stiffened elastic constant is

\[c^{p}_{KL}=c^{E}_{KL} + \frac{[e_{Kj} l_j][l_i e_{iL}]}{l_i \epsilon^{S}_{ij} l_j}.\]

We can explore this for different crystal cuts with piezo effects both disabled and active.

mat_LiNbO3 = materials.make_material("LiNbO3_1973_Auld")
#mat_LiNbO3 = materials.make_material("LiNbO3aniso_2021_Steel")

mat_LiNbO3_x = mat_LiNbO3.copy()
mat_LiNbO3_x.disable_piezoelectric_effects()
mat_LiNbO3_x.set_orientation('x-cut')

print(mat_LiNbO3_x.elastic_properties())

Elastic properties of material LiNbO3_2023_Rodrigues
  Density:        4650.000 kg/m^3
  Ref. index:     2.2100+0.0000j
  Crystal class:  Trigonal
  Crystal group:  no sym
  Stiffness c_IJ:
 Stiffness c_IJ, unit: GPa.
   Voigt 4-tensor:
    [[198.8300  54.6400  68.2300  -0.0000   7.8300   0.0000]
     [ 54.6400 198.8300  68.2300   0.0000  -7.8300  -0.0000]
     [ 68.2300  68.2300 235.7100   0.0000   0.0000   0.0000]
     [ -0.0000   0.0000   0.0000  59.8600   0.0000  -7.8300]
     [  7.8300  -7.8300   0.0000   0.0000  59.8600  -0.0000]
     [  0.0000  -0.0000   0.0000  -7.8300  -0.0000  72.0950]]

  Wave mode 1: v_p=7.1197 km/s,  |v_g|=7.1197 km/s,  u_j=[ 0.0000  0.0000  1.0000],  v_g=[ 0.0000  0.0000  7.1197] km/s
  Wave mode 2: v_p=3.5879 km/s,  |v_g|=3.6185 km/s,  u_j=[ 1.0000  0.0000  0.0000],  v_g=[ 0.4693 -0.0000  3.5879] km/s
  Wave mode 3: v_p=3.5879 km/s,  |v_g|=3.6185 km/s,  u_j=[ 0.0000  1.0000  0.0000],  v_g=[-0.4693  0.0000  3.5879] km/s

7.3.1.1. Piezo effects disabled

Starting with piezo effects disabled, the slowness curve (top right) in the first plot below corresponds to Fig. 8.9 in Auld volume 1.

imgoffx = mat_LiNbO3_x.plot_bulk_dispersion_all("tmp_xcut_piezo_off", show_poln=False, flip_x=True)
img_pair(imgoffx, "images/auld_fig8_9.png",  'X-cut, no piezo - NumBAT', 'X-cut, no piezo - Auld Fig. 8.9', 375, 300)

X-cut, no piezo - NumBAT

X-cut, no piezo - Auld Fig. 8.9

To make this comparison clearer, we can plot just the slowness curve (inverse phase velocity) in the top right. The flip_x=True flag reverses the sense of the \(x\)-axis. This compensates for NumBAT plotting dispersion in the \(q_x\)-\(q_z\) plane, while Auld uses the \(q_x\)-\(q_y\) plane.

imgoff2 = mat_LiNbO3_x.plot_bulk_dispersion_ivp("tmp_xcut_piezo_offb", show_poln=False, flip_x=True)
img_pair(imgoff2, "images/auld_fig8_9.png",  'X-cut, piezo off - NumBAT', 'X cut, piezo off - Auld Fig. 8.9', 420, 300)

X-cut, piezo off - NumBAT

X cut, piezo off - Auld Fig. 8.9

7.3.1.2. Piezo effects enabled

mat_LiNbO3_x.enable_piezoelectric_effects()
imgonx = mat_LiNbO3_x.plot_bulk_dispersion_ivp("tmp_xcut_piezo_on", show_poln=True, flip_x=True)
img_pair(imgonx, "images/auld_fig8_11.png",  'X-cut, piezo on - NumBAT', 'X-cut, piezo on - Auld Fig. 8.11', 400, 300)

X-cut, piezo on - NumBAT

X-cut, piezo on - Auld Fig. 8.11

The slowness curve in the second plot with piezo effects enabled corresponds to Fig. 8.11 in Auld volume 1.

mat_LiNbO3_z = mat_LiNbO3.copy()
mat_LiNbO3_z.set_orientation('z-cut')
#mat_LiNbO3_z.disable_piezoelectric_effects()

imgonz = mat_LiNbO3_z.plot_bulk_dispersion_ivp("tmp_zcut_piezo_on", show_poln=True, flip_x=True)
img_pair(imgonz, "images/auld_fig8_12.png",  'Z-cut, piezo on - NumBAT', 'Z-cut, piezo on - Auld Fig. 8.12', 400, 300)

Z-cut, piezo on - NumBAT

Z-cut, piezo on - Auld Fig. 8.12

This should give Auld Fig. 8.11.

mat_LiNbO3_negx.enable_piezoelectric_effects()
imgon = mat_LiNbO3_negx.plot_bulk_dispersion("tmp_negxcut_piezo_on", show_poln=False)
display(Image(imgon, width=500))

This should give Auld Fig. 8.12.

mat_LiNbO3_z = mat_LiNbO3.copy()
mat_LiNbO3_z.set_orientation('z-cut')

mat_LiNbO3_z.enable_piezoelectric_effects()
imgon = mat_LiNbO3_z.plot_bulk_dispersion("tmp_negxcut_piezo_on", show_poln=False)
display(Image(imgon, width=500))

7.4. Other material comparisons

Cubic material: GaAs

mat_GaAs = materials.make_material("GaAs_1970_Auld")

The slowness plot (top-right) below reproduces Fig. 7.2 of Auld vol. 1, while the ray surface plot (lower-left) reproduces Fig. 7.9.

mat_GaAs.disable_piezoelectric_effects()

fimg = mat_GaAs.plot_bulk_dispersion("tmp_defcut", show_poln=False)
display(Image(fimg, width=500))

Now let’s reproduce Fig. 7.3. To preserve the \(z\) axis, rotate \([1,-1,0]\) into the \(x\) axis and \([1,-1,1]\) into the horizontal plane, we need a \(\pi/4\) rotation around \(z\):

mat_rot = mat_GaAs.copy()
mat_rot.rotate([0,0,1], np.pi/4)
fimg = mat_rot.make_crystal_axes_plot('ttx')
display(Image(fimg, width=300))

fimg = mat_rot.plot_bulk_dispersion("tmp_defcut", show_poln=False)
display(Image(fimg, width=800))

Hexagonal material: Cadmium sulfide

mat_CdS = materials.make_material("CdS_1973_Auld")
mat_CdS.set_orientation('z-cut')
#print(mat_CdS.full_str())

We rotate to the \(z\)-cut orientation to get a dispersion cut through the horizontal symmetry plane.

The slowness curve (top-right) then reproduces Auld Fig. 7.6.

mat_CdS.set_orientation('z-cut')
fimg = mat_CdS.plot_bulk_dispersion("tmp_defcut", show_poln=False)
display(Image(fimg, width=600))