# 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)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.
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)