Second harmonic generation with LBO crystal

LBO is a biaxial crystal. That is it has three distinct principle refractive indices, \(n_x\), \(n_y\) and \(n_z\) with \(n_x<n_y<n_z\). The common practice is to confine the propagation vector \(\vec{k}\) in one of the three principle planes; xy, yz or zx. Under these configurations, LBO works exactly like a uniaxial crystal.

# Calculate RI
def getRI(wavelength, temperature):
    """
    wavelength : um
    temperature : degC
    Returns:
        n_x, n_y, n_z
    """

    delta_T = temperature - 20

    # Sellmeier equations
    n_x = np.sqrt(
        2.4542
        + 0.01125 / (wavelength**2 - 0.01135)
        - 0.01388 * wavelength**2
    )

    n_y = np.sqrt(
        2.5390
        + 0.01277 / (wavelength**2 - 0.01189)
        - 0.01849 * wavelength**2
        + 4.3025e-5 * wavelength**4
        - 2.9131e-5 * wavelength**6
    )

    n_z = np.sqrt(
        2.5865
        + 0.01310 / (wavelength**2 - 0.01223)
        - 0.01862 * wavelength**2
        + 4.5778e-5 * wavelength**4
        - 3.2526e-5 * wavelength**6
    )

    # Temperature corrections
    delta_nx = (
        (-3.76 * wavelength + 2.30)
        * 1e-6
        * (delta_T + 29.13e-3 * delta_T**2)
    )

    delta_ny = (
        (6.01 * wavelength - 19.40)
        * 1e-6
        * (delta_T - 32.89e-4 * delta_T**2)
    )

    delta_nz = (
        (1.50 * wavelength - 9.70)
        * 1e-6
        * (delta_T - 74.49e-4 * delta_T**2)
    )

    n_x += delta_nx
    n_y += delta_ny
    n_z += delta_nz

    return np.real(n_x), np.real(n_y), np.real(n_z)

Type I phase matching

For \(\theta\)= 90 \(^\circ\), the wavevector is confined to the xy plane. The polarization direction perpendicular to the plane would be along z axis, and hence the ordinary refractive index would be \(n_{z}\) as this will remain unchanged with the orientation of the propagation vector in the xy plane. Further \(\phi\)=0 \(^\circ\) confines it in the x direction. The extraordinary refractive index would be \(n_y\).

For this specific crystal orientation, the crystal is essentially a negative uniaxial crystal, where \(n_{e}<n_{o}\). For type I phase matching the following scheme is employed. Two extraordinary fundamental beams (refractive index \(n_y\) are used to one frequency doubled beam with ordinary refractive index (refractive index \(n_z\)).

\[ n_z+n_z \rightarrow n_y\]

# Type I phase mismatch function 
def phase_mismatch_typeI(wavelength, temperature,phi):
    """
    Type-I SHG:
    z + z -> y

    Output in meter inverse
    """
    lambda_shg = wavelength / 2
    phi= np.deg2rad(phi)
    n_x_shg, n_y_shg, _ = getRI(lambda_shg, temperature)
    _, _, n_z_fundamental = getRI(wavelength, temperature)

    # Ordinary refractive index
    no=n_z_fundamental
    ne=1/(np.sqrt((np.cos(phi)/n_y_shg)**2+(np.sin(phi)/n_x_shg)**2))
    return 1e6*(4*np.pi/wavelength)*(no - ne)

Type II phase matching

For type II-phase matching, the phase matching is performed in XZ plane, \(\phi\)= 0 \(^\circ\) and the ordinary refractive index would be \(n_y\). The crystal is again cut with \(\theta\)= 0 \(^\circ\). That would confine the propagation vector to z direction. The extraordinary refractive index would be \(n_x\).

\[ n_{x}+n_{y} \rightarrow n_{x}\]

# Type II phase mismatch function
def phase_mismatch_typeII(wavelength, temperature,theta):
    """
    Type-II SHG:
    x + y -> x

    Output in m inverse 
    """
    theta= np.deg2rad(theta)
    lambda_shg = wavelength / 2

    n_x_shg, _, n_z_shg = getRI(lambda_shg, temperature)
    n_x_fundamental, n_y_fundamental, n_z_fundamental = getRI(wavelength, temperature)

    no=n_y_fundamental
    ne_fundamental=1/(np.sqrt((np.cos(theta)/n_x_fundamental)**2+(np.sin(theta)/n_z_fundamental)**2))
    ne_shg=1/(np.sqrt((np.cos(theta)/n_x_shg)**2+(np.sin(theta)/n_z_shg)**2))
    return 1e6*(-2 * ne_shg + (no + ne_fundamental))*(2*np.pi/wavelength)

Acceptance bandwidth calculation

# PM temperature calculation function

def PM_temp(wavelength,PM_type):
    temperatures = []

    if PM_type=="Type-I":
        for lam in wavelength:
            sol = root_scalar(
                lambda T: phase_mismatch_typeI(lam, T,phi=0),
                bracket=[-20, 300],
                method="brentq",)

            temperatures.append(sol.root)
    if PM_type=="Type-II":
        for lam in wavelength:
            sol = root_scalar(
                lambda T: phase_mismatch_typeII(lam, T,theta=0),
                bracket=[-20, 50],
                method="brentq",)

            temperatures.append(sol.root)
    temperatures = np.array(temperatures)
    return temperatures
        
# Angular acceptance bandwidth calculation for type-I
wavelengths=np.linspace(1.06,1.198,50)
temperatures=PM_temp(wavelength=wavelengths,PM_type="Type-I")

L=25e-3         # LBO crystal length in meters
acc_bndwd=[]   # Phase matching bandwidth in um
PM_phi=[]

# for i in range(len(wavelength1)):
for i in range(len(wavelengths)):
    sol = root_scalar(lambda p: phase_mismatch_typeI(wavelengths[i], 25, phi=p),bracket=[0, 90], method="brentq")
    PM_phi_val = sol.root
    phi_scan = np.linspace(PM_phi_val-0.8, PM_phi_val+0.8, 1000)

    dk = np.array([phase_mismatch_typeI(wavelengths[i], 25,phi=t) for t in phi_scan])

    

    
    eff = (np.sinc(dk*L/(2*np.pi)))**2
    half_value = np.max(eff)*0.5
    spline = UnivariateSpline(phi_scan, eff-half_value, s=0)
    roots = spline.roots()  # Returns points where the spline crosses 0
    
    acc_bndwd_val = roots[1] - roots[0]
    
    acc_bndwd.append(acc_bndwd_val)
    PM_phi.append(PM_phi_val)
    
acc_bndwd = np.array(acc_bndwd)
acc_bndwd_mrad=np.deg2rad(acc_bndwd)*1e3
PM_phi=np.array(PM_phi)

# Phase matching temperature calculation
# Type-I
wavelength1 = np.linspace(1.06, 1.6, 50)
temperatures1 = PM_temp(wavelength=wavelength1,PM_type="Type-I")


# Type-II
wavelength2 = np.linspace(1.15, 1.45, 50)

temperatures2 =PM_temp(wavelength=wavelength2,PM_type="Type-II")


L=25e-3         # LBO crystal length in meters
linewidth1=[]   # Phase matching bandwidth in um
wavelength1=np.linspace(1.060,1.2,100)  # Wavelengths in um
temperatures1=PM_temp(wavelength=wavelength1,PM_type="Type-I")
# for i in range(len(wavelength1)):
for i in range(len(wavelength1)):
    lambda1=wavelength1[i]
    T=temperatures1[i]
    # print(lambda1)
    # print(T)
    lam_scan = np.linspace(lambda1-0.005, lambda1+0.005, 1000)

    dk = np.array([phase_mismatch_typeI(l, T,phi=0) for l in lam_scan])

    eff = (np.sinc(dk*L/(2*np.pi)))**2
    half_value = np.max(eff)*0.5
    spline = UnivariateSpline(lam_scan, eff-half_value, s=0)
    roots = spline.roots()  # Returns points where the spline crosses 0
    linewidth_val = roots[1] - roots[0]
    # print(len(roots))
    linewidth1.append(linewidth_val)
    
linewidth1 = np.array(linewidth1)

linewidth2=[]
wavelength2=np.linspace(1.2,1.26,100)
temperatures2=PM_temp(wavelength=wavelength2,PM_type="Type-II")
# for i in range(len(wavelength1)):
for i in range(len(wavelength2)):
    lambda1=wavelength2[i]
    T=temperatures2[i]
    # print(lambda1)
    # print(T)
    lam_scan = np.linspace(lambda1-0.05, lambda1+0.05, 1000)

    dk = np.array([phase_mismatch_typeII(l, T,theta=0) for l in lam_scan])

    eff = (np.sinc(dk*L/(2*np.pi)))**2
    half_value = np.max(eff)*0.5
    spline = UnivariateSpline(lam_scan, eff-half_value, s=0)
    roots = spline.roots()  # Returns points where the spline crosses 0
    linewidth_val = roots[1] - roots[0]
    linewidth2.append(linewidth_val)
    # plt.plot(eff)

    
linewidth2 = np.array(linewidth2)

# plt.ylim(0,10)