Problem with phonon simulation

[LaureTailpied]

Hello, Thank you very much for the development of the package, it is very useful. We are trying to use it to simulate the response of phonon polaritons (for instance for h-BN), but it gives strange results (see picture enclosed) Do you have an idea on how we can improve the simulation, or if it is a result of the coding itself?

Thank you very much for your help, Best regards

Laure Tailpied

[TomVincentUK]

Hi Laure,

Thanks for using our package! It's a bit tricky to know what the issue is just from the attached picture, but I'd caution that simulating polaritons can be tricky as they involve strong light-matter interactions. Also, for hBN there's a strong anisotropy to the permittivity which we can't account for yet in snompy. If you're happy to share your script to generate the figure, and to give me some indication of how you'd expect the figure to look, I can take a look to see if there's anything obviously at fault, but I can't promise that I'll be able to see anything.

Cheers, Tom

On Tue, 2 Jul 2024, 9:54 am LaureTailpied, @.***> wrote:

Hello, Thank you very much for the development of the package, it is very useful. We are trying to use it to simulate the response of phonon polaritons (for instance for h-BN), but it gives strange results (see picture enclosed) 2024_07_Pb_Phonons.png (view on web) https://github.com/TomVincentUK/snompy/assets/174425694/a8b2e7f7-a263-409b-b007-79915170bfb9 Do you have an idea on how we can improve the simulation, or if it is a result of the coding itself?

Thank you very much for your help, Best regards

Laure Tailpied

— Reply to this email directly, view it on GitHub https://github.com/TomVincentUK/snompy/issues/41, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANNTHFVFVEKPUECMRWQGPFDZKJTELAVCNFSM6AAAAABKHCVOGKVHI2DSMVQWIX3LMV43ASLTON2WKOZSGM4DKNRVGMYDENQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[LaureTailpied]

Hi Tom, Thank you very much for your answer. Attached is my code for the simulation of GaSb, the phase at least should be more continuous with a single peak. `# -- coding: utf-8 -- """ Created on Thu Jun 13 15:08:51 2024

@author: bfix modified by ltailpie """

import matplotlib.pyplot as plt import numpy as np import snompy

c = 31e8 #m/s speed of light pi = np.pi eps_0 = 8.85418782e-12 m0 = 9.109e-31 #masse electron e = 1.61e-19 #Charge electron

%% Set some experimental parameters

A_tip = 65e-9 # AFM tip tapping amplitude r_tip = 10e-9 # AFM tip radius of curvature L_tip = 200e-9 # Semi-major axis length of ellipsoid tip model n = 3 # Harmonic for demodulation theta_in = np.deg2rad(60) # Light angle of incidence c_r = 0.3 # Experimental weighting factor nu_vac = np.linspace(500, 50 -00, 4096) * 1e2 # Vacuum wavenumber method = "Q_ave" # The FDM method to use

%%Fonctions diélectriques

Semi-infinite superstrate and substrate

eps_air = 1.0 eps_Si = 11.7 # Si permitivitty in the mid-infrared

def eps_GaSb(nu_vac,model='Kukharskii') :

Kukharskii mobility of doped semiconductors extracted from Fourier-transform infrared ellipsometry spectra

# Stefan Zollner ; Pablo P. Paradis; Farzin Abadizaman; Nuwanjula S. Samarasingha
# JVST B (2019)
# Et SI
# https://doi.org/10.1116/1.5081055
lambda_ = 1/nu_vac
E_eV = 1.24/(lambda_*1e6)
# dopage p 2.5e18cm-3
# screened plasma energy
eps_infini = 14
E_LP = 0.0251
Gamma_LP = 0.001
Gamma_K = 0.0064
E_UP = 0.032
Gamma_UP = 0.0117
E_TO = 0.02778
Gamma_TO = 0.0003
eps_pGaSb = eps_infini * ((E_LP**2-E_eV**2-1j*Gamma_LP*E_eV)/(-E_eV*(E_eV+1j*Gamma_K))) * ((E_UP**2-E_eV**2-1j*Gamma_UP*E_eV)/(E_TO**2-E_eV**2-1j*Gamma_TO*E_eV));
return(eps_pGaSb)

eps_GaSb = eps_GaSb(nu_vac)

%%Création échantillon

sample_GaSb = snompy.bulk_sample( eps_GaSb) sample_Si = snompy.bulk_sample(eps_sub=11.7, eps_env=eps_air, nu_vac=nu_vac)

%% Measurement

alpha_eff_GaSb = snompy.fdm.eff_pol_n( sample=sample_GaSb, A_tip=A_tip, n=n, r_tip=r_tip, L_tip=L_tip, method=method ) r_coef_GaSb = sample_GaSb.refl_coef(theta_in=theta_in) sigma_GaSb = (1 + c_r * r_coef_GaSb) * 2 alpha_eff_GaSb

Si reference

alpha_eff_Si = snompy.fdm.eff_pol_n( sample=sample_Si, A_tip=A_tip, n=n, r_tip=r_tip, L_tip=L_tip, method=method ) r_coef_Si = sample_Si.refl_coef(theta_in=theta_in) sigma_Si = (1 + c_r * r_coef_Si) * 2 alpha_eff_Si

Normalised complex scattering

eta_GaSb = sigma_GaSb / sigma_Si angle_GaSb = np.angle(sigma_GaSb)-np.angle(sigma_Si)

%% Plot output GaSb

fig, axes = plt.subplots(nrows=2, sharex=True)

For neater plotting

nu_per_cm = nu_vac * 1e-2

axes[0].plot(nu_per_cm, np.abs(eta_GaSb), c='C0',label='Charge average method') axes[1].plot(nu_per_cm,(angle_GaSb), c='C1',label='Charge average method')

axes[0].setylabel(r"$s{" f"{n}" r"}$ / a.u.") axes[1].set( xlabel=r"$\nu$ / cm$^{-1}$", ylabel=r"$\phi_{" f"{n}" r"}$ / radians", xlim=(nu_per_cm.max(), nu_per_cm.min()), ) axes[0].legend() axes[1].legend() fig.tight_layout() plt.show(block=False)`

And here is the code for the h-BN, which should be more linear + have a single peak in phase at ~1300 cm-1 in place of the inversion

""" Created on Thu Jun 13 15:08:51 2024

@author: bfix """

import matplotlib.pyplot as plt import numpy as np import snompy

c = 31e8 #m/s speed of light pi = np.pi eps_0 = 8.85418782e-12 m0 = 9.109e-31 #masse electron e = 1.61e-19 #Charge electron

Set some experimental parameters

A_tip = 80e-9 # AFM tip tapping amplitude r_tip = 10e-9 # AFM tip radius of curvature L_tip = 600e-9 # Semi-major axis length of ellipsoid tip model n = 3 # Harmonic for demodulation theta_in = np.deg2rad(60) # Light angle of incidence c_r =0.3 # Experimental weighting factor nu_vac = np.linspace(700, 1800, 1024) * 1e2 # Vacuum wavenumber method = "Q_ave" # The FDM method to use

Semi-infinite superstrate and substrate

eps_air = 1.0 eps_Si = 11.7 # Si permitivitty in the mid-infrared

def eps_hBN(nu_vac): omega = 2pic*nu_vac

perpendicular phonon polaritons only

epsiinfPP=4.87
epsistatPP=6.71
SPP=epsistatPP-epsiinfPP
wtoPP=1369.8*2*pi*c*100
wtoPP=1363*2*pi*c*100
gamPP=5*2*pi*c*100 # damping, taken as 3% of wTO
epsiPP=epsiinfPP+(SPP*wtoPP**2)/(wtoPP**2-omega**2-1j*gamPP*omega)
epsi1PP = np.real(epsiPP)
epsi2PP = np.imag(epsiPP)
# parallel phonon polaritons only
epsiinfP=2.95
epsistatP=3.57
SP=epsistatP-epsiinfP;
wtoP=800*2*pi*c*100
gamP=0.2*2*pi*c*100 #damping, taken as 3% of wTO
epsiP=epsiinfP+(SP*wtoP**2)/(wtoP**2-omega**2-1j*gamP*omega)
epsi1P=np.real(epsiP)
epsi2P=np.imag(epsiP)
return ((epsi1P + 1j*epsi2P) + 2*(epsi1PP + 1j*epsi2PP))/3

def eps_SiO2(nu_vac): omega = 2picnu_vac epsiinf3=1.97 wto3=10752pic100 wlo3=12502pic100 epsistat3=epsiinf3((wlo32)/(wto32)) S3=epsistat3-epsiinf3 gam3=692pic100 # damping, epsi3 =epsiinf3+S3*wto32/(wto32-omega*2-1jgam3*omega) return epsi3

eps_hbn = eps_hBN(nu_vac) eps_SiO2 = eps_SiO2(nu_vac) sample_BN = snompy.Sample( eps_stack=(eps_air, eps_hbn, eps_SiO2), t_stack=(1*1e-9,), nu_vac=nu_vac)

Model of Au dielectric function from ref [2] below

sample_Si = snompy.bulk_sample(eps_sub=11.7, eps_env=eps_air, nu_vac=nu_vac)

Measurement

alpha_eff_BN = snompy.fdm.eff_pol_n( sample=sample_BN, A_tip=A_tip, n=n, r_tip=r_tip, L_tip=L_tip, method=method ) r_coef_BN = sample_BN.refl_coef(theta_in=theta_in) sigma_BN = (1 + c_r * r_coef_BN) * 2 alpha_eff_BN

Si reference

alpha_eff_Si = snompy.fdm.eff_pol_n( sample=sample_Si, A_tip=A_tip, n=n, r_tip=r_tip, L_tip=L_tip, method=method ) r_coef_Si = sample_Si.refl_coef(theta_in=theta_in) sigma_Si = (1 + c_r * r_coef_Si) * 2 alpha_eff_Si

Normalised complex scattering

eta_BN = sigma_BN / sigma_Si angle_BN = np.unwrap(np.angle(alpha_eff_BN)-np.angle(alpha_eff_Si))

%% Plot output BN

fig, axes = plt.subplots(nrows=2, sharex=True)

For neater plotting

nu_per_cm = nu_vac * 1e-2

axes[0].plot(nu_per_cm, np.abs(eta_BN), c='C0') axes[1].plot(nu_per_cm,np.angle(eta_BN), c='C1')

axes[0].setylabel(r"$s{" f"{n}" r"}$ / a.u.") axes[1].set( xlabel=r"$\nu$ / cm$^{-1}$", ylabel=r"$\phi_{" f"{n}" r"}$ / radians", xlim=(nu_per_cm.max(), nu_per_cm.min()), ) fig.tight_layout()

plt.show(block=False)

Thank you very much for your help ! Laure

[TomVincentUK]

Hey Laure,

I think the phase discontinuity you see is just an unwrapping issue (because phi + n2pi is the same as phi). You can solve that by replacing your (angle_GaSb) from the plotting line with np.unwrap(angle_GaSb). (Also you don't need to calculate the angle by subtracting the reference phase from the sample phase as you did here, you can just use np.angle(eta_GaSb) for the same thing).

I'll repeat my cautions about strong light-matter interactions and the anisotropy of hBN, which might lead to apparently strange results.

I'm going to mark this Issue as closed, because I think this is more of a question about how to use these models, rather than a problem with the code. But if you've got more questions on this topic, please feel free to start a new discussion under the Discussions tab of this repository.

Cheers, Tom