دفترچهٔ فصل ۳: کتابخانهٔ معتبر سلول دی‌الکتریک

شناسه: ms-07-unit-cell. محاسبهٔ ماکسول تمام‌موج حوزهٔ فرکانس برای توری ریج دی‌الکتریک یک‌بعدی، بی‌اتلاف و غیرمغناطیسی با تابش TE، یعنی $E_y,k_y=0$ ، با ماژول اصیل simulations/unit_cell_te.py، NumPy، SciPy و Matplotlib اجرا می‌شود. بستهٔ منبع یا کلون مخزن، ماژول را حفظ می‌کند.

پیش‌بینی: آیا عرض تنها، همهٔ فازها را با عبور زیاد می‌دهد؟ پاسخ فقط برای همین هندسهٔ ایده‌آل است؛ ستون استوانه‌ای، پاشندگی، ساخت، آرایهٔ متناهی و قطبش دلخواه نیست.

قرارداد $e^{-i\omega t}$ ، $e^{i\mathbf{k}\cdot\mathbf{r}}$ ، تابش از $z<0$ ، ریج مرکزی $|x|<w/2,0<z<h$ و واحد طول میکرومتر است. $E_y$ فرودی در $z=0$ واحد؛ مرجع $r_m$ بالا و $t_m$ در $z=h$ است. انتشار عادی نتایج ذخیره‌شده را نمایش می‌دهد و حل‌گر اجرا نمی‌کند.

from pathlib import Path
import sys, json, csv, hashlib, platform, html
from dataclasses import replace, asdict
import numpy as np
import scipy
import matplotlib
import matplotlib.pyplot as plt
from IPython.display import display, HTML

ROOT = next((p for p in [Path.cwd(), *list(Path.cwd().parents)[:3]]
             if (p/'simulations/unit_cell_te.py').is_file()), None)
if ROOT is None:
    raise FileNotFoundError('Extract the Chapter 3 source ZIP, or run from the repository notebooks folder; simulations/unit_cell_te.py is required.')
sys.path.insert(0, str(ROOT/'simulations'))
from unit_cell_te import Cell, solve_cell, film_reference, deembed_te, phase_error_deg, sample_ey, SOLVER_VERSION

# EDIT THESE PARAMETERS. Changing them defines a new experiment.
REFERENCE = Cell(period_um=0.65, width_um=0.325, height_um=0.60,
                 wavelength_um=1.55, n_ridge=2.4, n_gap=1.0,
                 n_incident=1.0, n_substrate=1.5, angle_deg=0.0)
WIDTHS_UM = np.round(np.linspace(0.08, 0.57, 50), 8)
M_LIBRARY, M_CHECK = 24, 32
POWER_TOL, COMPLEX_TOL, PHASE_TOL_DEG = 5e-4, 5e-4, 0.05
MIN_T_FOR_SELECTION, TARGET_TOL_DEG = 0.8, 5.0
ENABLE_WIDGETS = False
DATA = ROOT/'data/chapter-03'
DATA.mkdir(parents=True, exist_ok=True)

def table(headers, rows):
    def escaped(value):
        return html.escape(f'{value:.6g}' if isinstance(value,(float,np.floating)) else str(value))
    display(HTML('<table><thead><tr>'+''.join('<th>'+escaped(x)+'</th>' for x in headers)+
                 '</tr></thead><tbody>'+''.join('<tr>'+''.join('<td>'+escaped(x)+'</td>' for x in row)+'</tr>' for row in rows)+'</tbody></table>'))

plt.rcParams.update({'figure.dpi': 125, 'font.size': 10, 'axes.spines.top': False, 'axes.spines.right': False})
print('Python', platform.python_version(), '| NumPy', np.__version__, '| SciPy', scipy.__version__, '| Matplotlib', matplotlib.__version__)
print('Solver:', SOLVER_VERSION, '| harmonics:', 2*M_LIBRARY+1, 'checked against', 2*M_CHECK+1)
table(['Parameter', 'Value'], asdict(REFERENCE).items())

Python 3.12.1 | NumPy 2.4.6 | SciPy 1.18.1 | Matplotlib 3.11.1
Solver: te-ridge-1.0 | harmonics: 49 checked against 65

۱. پاسخ معلوم مستقل را بررسی کنید¶

در $h=0$ ، مرز تخت با $r=(q_i-q_t)/(q_i+q_t)$ و $t=2q_i/(q_i+q_t)$ داریم. در $w=0,w=\Lambda$ لایه یکنواخت است. فرمول فابری‌پروی مستقل دامنهٔ مختلط می‌دهد، نه فقط $R+T$ . دو حد لایه را در سه زاویه و سه ارتفاع کنترل کنید. این‌ها مرجع و علامت را می‌سنجند؛ ریج نقش‌دار هنوز همگرایی مرتبه و مقایسهٔ جدا می‌خواهد.

benchmark_rows=[]
max_complex_error=0.0
for angle in [0.,20.,45.]:
    for height in [0.,0.17,0.60]:
        for width, index in [(0.,REFERENCE.n_gap),(REFERENCE.period_um,REFERENCE.n_ridge)]:
            cell=replace(REFERENCE,angle_deg=angle,height_um=height,width_um=width)
            result=solve_cell(cell,8)
            ra,ta=film_reference(cell,index)
            error=max(abs(result['r'][8]-ra),abs(result['t'][8]-ta))
            max_complex_error=max(max_complex_error,error)
            assert error < 1e-10
            assert result['balance_error'] < 1e-10
        benchmark_rows.append([angle,height,error,result['balance_error']])
flat=solve_cell(replace(REFERENCE,height_um=0),8)
assert abs(flat['r'][8]+0.2)<1e-10 and abs(flat['t'][8]-0.8)<1e-10
print('18 uniform-film/interface checks passed; maximum complex-amplitude error:',max_complex_error)
print('Normal flat interface: r=',flat['r'][8], 't=',flat['t'][8], 'R=',flat['R_total'],'T=',flat['T_total'])
table(['Angle (deg)','Height (um)','Full-width film complex error','Power-balance error'],benchmark_rows)

18 uniform-film/interface checks passed; maximum complex-amplitude error: 8.865462800015876e-15
Normal flat interface: r= (-0.20000000000000007+1.2004609669086259e-17j) t= (0.7999999999999999+1.2004609669086259e-17j) R= 0.04000000000000003 T= 0.9599999999999999

۲. باقیماندهٔ کوچک توان، همگرایی نیست¶

مسئلهٔ ویژهٔ TE برابر $[\epsilon_{m-n}-\alpha_m^2\delta_{mn}]W=WQ^2$ با $\alpha_m=n_i\sin\theta_i+m\lambda_0/\Lambda$ است. مستطیل ضرایب دقیق $\epsilon_\ell=n_g^2\delta_{\ell0}+(n_r^2-n_g^2)f\,\mathrm{sinc}(\ell f)$ با $f=w/\Lambda$ و $\mathrm{sinc}(u)=\sin(\pi u)/(\pi u)$ دارد. میدان و مشتق عمودی روی دو وجه تطبیق می‌یابند.

در عرض ۰٫۳۲۵ میکرومتر، نیم‌مرتبه‌های $M$ را با $M=48$ مقایسه کنید. $M=0$ می‌تواند باقیماندهٔ کوچک و عبور مختلط نادقیق داشته باشد. کتابخانه همهٔ عرض‌ها را در $M=24,32$ کنترل می‌کند.

reference_high=solve_cell(REFERENCE,48)
convergence=[]
for M in [0,1,2,4,8,12,16,24,32]:
    result=solve_cell(REFERENCE,M)
    convergence.append([M,2*M+1,result['T_total'],
                        abs(result['t'][M]-reference_high['t'][48]),
                        abs(phase_error_deg(result['t'][M],reference_high['t'][48])),
                        result['balance_error']])
table(['M','Harmonics','T total','Complex t error vs M=48','Phase error (deg)','Power residual'],convergence)
assert convergence[-1][3] < COMPLEX_TOL

۳. میدان سازندهٔ پاسخ را ببینید¶

ضرایب مدی میدان فرودی، بازتابی، داخلی و عبوری را بازسازی می‌کنند. خطوط سیاه ریج تناوبی و زیرلایه‌اند. $|E_y/E_{\rm inc}|^2$ مربع دامنه است، نه توان عبوری یا چگالی انرژی. هارمونیک بسته با توان دورمیدان صفر، در نزدیکی ریج سهم دارد.

result=solve_cell(REFERENCE,M_LIBRARY)
x=np.linspace(-1.5*REFERENCE.period_um,1.5*REFERENCE.period_um,281)
z=np.linspace(-0.45,REFERENCE.height_um+0.65,241)
field=sample_ey(REFERENCE,result,x,z)
fig,axes=plt.subplots(1,2,figsize=(11,4.5),layout='constrained')
for ax,values,title,cmap in zip(axes,[field.real,abs(field)**2],['Real Ey / Einc','Squared field amplitude |Ey / Einc|²'],['RdBu_r','magma']):
    image=ax.pcolormesh(x,z,values,shading='auto',cmap=cmap)
    fig.colorbar(image,ax=ax,shrink=0.85)
    for center in [-REFERENCE.period_um,0,REFERENCE.period_um]:
        ax.add_patch(plt.Rectangle((center-REFERENCE.width_um/2,0),REFERENCE.width_um,REFERENCE.height_um,fill=False,lw=1.1,color='black'))
    ax.axhline(REFERENCE.height_um,color='black',lw=1)
    ax.set(xlabel='x (um)',ylabel='z (um)',title=title)
    ax.invert_yaxis()
fig.suptitle('TE ridge unit cell: full-wave Fourier reconstruction; incidence travels downward')
plt.show()

۴. کتابخانه و شاهد آن را حفظ کنید¶

در این دوره و طول موج فقط $m=0$ در دو بیرون باز است. در قرارداد TE، $R_m=\mathrm{Re}(q_{i,m})|r_m|^2/q_{i,0}$ و $T_m=\mathrm{Re}(q_{t,m})|t_m|^2/q_{i,0}$ ؛ ریشه‌ها حقیقی مثبت یا موهومی مثبت‌اند.

هر ردیف حقیقی و موهومی، توان، باقیمانده و تغییر با $M$ را نگه می‌دارد. فراداده، هندسه، قرارداد، صفحات مرجع، هش، محیط، آستانه و دامنهٔ اعتبار دارد. این سلول JSON و CSV را بازمی‌سازد؛ فقط ردیف گذرنده پذیرفته است. مقایسهٔ قطع‌های متناهی کران خطای ریاضی تضمین‌شده نیست.

rows=[]
for width in WIDTHS_UM:
    cell=replace(REFERENCE,width_um=float(width))
    low,high=solve_cell(cell,M_LIBRARY),solve_cell(cell,M_CHECK)
    t,r=low['t'][M_LIBRARY],low['r'][M_LIBRARY]
    dt=abs(t-high['t'][M_CHECK])
    dp=max(abs(low['T_total']-high['T_total']),abs(low['R_total']-high['R_total']))
    dphi=abs(phase_error_deg(t,high['t'][M_CHECK]))
    accepted=bool(dt<COMPLEX_TOL and dp<POWER_TOL and dphi<PHASE_TOL_DEG and low['balance_error']<1e-8)
    rows.append(dict(width_um=float(width),r_real=float(r.real),r_imag=float(r.imag),
                     t_real=float(t.real),t_imag=float(t.imag),T=low['T_total'],R=low['R_total'],
                     phase_deg=float(np.rad2deg(np.angle(t))),balance_error=low['balance_error'],
                     solve_residual=low['solve_residual'],delta_t=float(dt),delta_power=float(dp),
                     delta_phase_deg=float(dphi),accepted=accepted))
assert all(row['accepted'] for row in rows)
t_values=np.array([complex(row['t_real'],row['t_imag']) for row in rows])
T_values=np.array([row['T'] for row in rows])
unwrapped=np.rad2deg(np.unwrap(np.angle(t_values)))
phase_span=float(np.ptp(unwrapped))
metadata=dict(schema='metasurface-unit-cell-library-v1',topic_id='ms-07-unit-cell',
              model='1D centered rectangular dielectric ridge, TE Ey, ky=0, single uniform-height layer',
              solver=SOLVER_VERSION,source_sha256=hashlib.sha256((ROOT/'simulations/unit_cell_te.py').read_bytes()).hexdigest(),
              parameters=asdict(REFERENCE),varying_parameter='width_um',units='micrometers and degrees',
              time_convention='exp(-i omega t)',spatial_convention='exp(i k.r)',
              phase_references={'r':'z=0','t':'z=height','incident_Ey':'one at z=0'},
              material_model='positive real constant indices; nonmagnetic, lossless, nondispersive',
              harmonics=2*M_LIBRARY+1,check_harmonics=2*M_CHECK+1,
              thresholds=dict(delta_complex_t=COMPLEX_TOL,delta_power=POWER_TOL,delta_phase_deg=PHASE_TOL_DEG,power_residual=1e-8),
              environment={'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__},
              validation_scope='uniform-film analytic limits, width-wise harmonic refinements; see validation.json for additional audits',
              unwrapped_phase_span_deg=phase_span,
              limitations=['No TM or conical incidence','No 2D lattice or finite array','No absorption or dispersion','No fabricated-device validation','No full 360-degree phase coverage for this width sweep'])
library=dict(metadata=metadata,rows=rows)
(DATA/'response-library.json').write_text(json.dumps(library,indent=2)+'\n',encoding='utf-8')
with (DATA/'response-library.csv').open('w',newline='',encoding='utf-8') as handle:
    writer=csv.DictWriter(handle,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows)
print(f'{len(rows)} / {len(rows)} widths passed. Continuous sampled phase span: {phase_span:.3f} degrees.')
print('Worst width-wise changes: complex t',max(row['delta_t'] for row in rows),'; power',max(row['delta_power'] for row in rows),'; phase (deg)',max(row['delta_phase_deg'] for row in rows))
table(['Width (um)','R','T','Phase (deg)','Delta complex t','Accepted'],
      [[row[k] for k in ['width_um','R','T','phase_deg','delta_t','accepted']] for row in rows[::7]])

50 / 50 widths passed. Continuous sampled phase span: 132.042 degrees.
Worst width-wise changes: complex t 6.5843872823157e-06 ; power 1.122202837988695e-06 ; phase (deg) 0.0004712439155618459

fig,axes=plt.subplots(1,3,figsize=(12,3.8),layout='constrained')
axes[0].plot(WIDTHS_UM,T_values,label='T');axes[0].plot(WIDTHS_UM,[r['R'] for r in rows],label='R')
axes[0].set(xlabel='Width (um)',ylabel='Incident-power fraction',ylim=(0,1.05),title='Power response');axes[0].legend()
axes[1].plot(WIDTHS_UM,[r['phase_deg'] for r in rows],'.',label='Wrapped')
axes[1].plot(WIDTHS_UM,unwrapped,label='Unwrapped')
axes[1].set(xlabel='Width (um)',ylabel='Transmission phase (deg)',title=f'Sampled continuous span: {phase_span:.1f}°');axes[1].legend()
axes[2].plot(t_values.real,t_values.imag,'o-',ms=3)
theta=np.linspace(0,2*np.pi,200);axes[2].plot(np.cos(theta),np.sin(theta),'--',color='0.6',lw=1)
axes[2].set(xlabel='Re(t)',ylabel='Im(t)',title='Complex field amplitude',aspect='equal',xlim=(-1.1,1.1),ylim=(-1.1,1.1))
for ax in axes: ax.grid(alpha=.2)
fig.suptitle('Width changes amplitude and phase together; it does not guarantee a complete phase library')
plt.show()

۵. همسایگی طیفی¶

کتابخانه یک طول موج و زاویه را پوشش می‌دهد. جاروب سه عرض در ۱٫۳۵ تا ۱٫۷۵ میکرومتر لزوم حفظ فاز و توان را نشان می‌دهد، ولی دامنهٔ معتبر صادرات را گسترش نمی‌دهد. کتابخانهٔ جدید در طول موج و زاویهٔ خودش کنترل هر ردیف می‌خواهد. ضریب ثابت، پاشندگی ماده را عمداً حذف می‌کند.

wavelengths=np.linspace(1.35,1.75,41)
fig,axes=plt.subplots(1,2,figsize=(10,3.7),layout='constrained')
for width in [.15,.325,.50]:
    results=[solve_cell(replace(REFERENCE,width_um=width,wavelength_um=float(w)),M_LIBRARY) for w in wavelengths]
    transmissions=np.array([v['t'][M_LIBRARY] for v in results])
    axes[0].plot(wavelengths,[v['T_total'] for v in results],label=f'w={width:g} um')
    axes[1].plot(wavelengths,np.rad2deg(np.unwrap(np.angle(transmissions))),label=f'w={width:g} um')
    assert max(v['balance_error'] for v in results)<1e-8
axes[0].set(ylabel='T',ylim=(0,1.05),title='Transmission spectrum')
axes[1].set(ylabel='Unwrapped phase (deg)',title='Phase spectrum at fixed reference planes')
for ax in axes: ax.set_xlabel('Free-space wavelength (um)');ax.grid(alpha=.2);ax.legend()
plt.show()

۶. انتخاب و گزارش صادقانهٔ درخواست ناممکن¶

برای $\phi_\star$ ، خطای دایروی $\Delta\phi=\arg(t e^{-i\phi_\star})$ را بگیرید. فقط پذیرفته‌ها با $T\geq0.8$ را نگه و $|\Delta\phi|$ را کمینه کنید. امکان‌پذیری به خطای درون $5^\circ$ نیاز دارد. کتابخانهٔ غیرخالی همیشه نزدیک‌ترین ردیف دارد؛ طرح رضایت‌بخش الزاماً ندارد. این انتخاب گسسته است، نه بهینه‌سازی کامل فراسطح.

eligible=[row for row in rows if row['accepted'] and row['T']>=MIN_T_FOR_SELECTION]
selection=[]
for target in [-160.,-100.,-50.,0.]:
    if not eligible:
        selection.append([target,'none','none','none','no eligible geometry']);continue
    candidate=min(eligible,key=lambda row:abs(phase_error_deg(complex(row['t_real'],row['t_imag']),np.exp(1j*np.deg2rad(target)))))
    error=abs(phase_error_deg(complex(candidate['t_real'],candidate['t_imag']),np.exp(1j*np.deg2rad(target))))
    selection.append([target,candidate['width_um'],candidate['T'],error,'feasible' if error<=TARGET_TOL_DEG else 'target unmet'])
table(['Target phase (deg)','Nearest eligible width (um)','T','Circular error (deg)','Decision'],selection)

۷. مانیتور، فاز را بدون تغییر قطعه عوض می‌کند¶

مانیتور ورودی در $z=-L_i$ و خروجی در $z=h+L_t$ ؛ نسبت به میدان فرودی مانیتور اول، $r_{\rm mon}=r e^{2ik_0q_iL_i}$ و $t_{\rm mon}=t e^{ik_0(q_iL_i+q_tL_t)}$ . سلول بعد جابه‌جایی را از دامنهٔ محاسبه‌شده می‌سازد و حذف می‌کند. آزمون جبری مرجع است؛ خود RCWA روی دو مرز و بدون مانیتور جابه‌جا تطبیق می‌دهد.

result=solve_cell(REFERENCE,M_LIBRARY)
r,t=result['r'][M_LIBRARY],result['t'][M_LIBRARY]
qi,qt=result['qi'][M_LIBRARY],result['qt'][M_LIBRARY]
k0=2*np.pi/REFERENCE.wavelength_um
monitor_rows=[]
for Li,Lt in [(0.,0.),(.1,.2),(.3,.4)]:
    rm=r*np.exp(2j*k0*qi*Li)
    tm=t*np.exp(1j*k0*(qi*Li+qt*Lt))
    rd,td=deembed_te(rm,tm,REFERENCE.wavelength_um,qi,qt,Li,Lt)
    assert max(abs(rd-r),abs(td-t))<1e-12
    monitor_rows.append([Li,Lt,np.rad2deg(np.angle(tm)),np.rad2deg(np.angle(td)),abs(td-t)])
table(['Li (um)','Lt (um)','Raw monitor t phase (deg)','Face-referenced phase (deg)','Complex error'],monitor_rows)

۸. کنترل اختیاری ژوپیتر¶

ENABLE_WIDGETS=True را در پارامتر تنظیم و دوباره اجرا کنید تا کتابخانهٔ محاسبه‌شده را ببینید. کنترل، ردیف موجود را انتخاب و مسئلهٔ تازه حل نمی‌کند. آزمایشگاه منتشرشده مرورگر مستقل بدون نصب پایتون نیز دارد.

if ENABLE_WIDGETS:
    import ipywidgets as widgets
    def inspect_width(index):
        row=rows[index]
        table(['Width (um)','T','Phase (deg)','Delta complex t'],[[row['width_um'],row['T'],row['phase_deg'],row['delta_t']]])
    display(widgets.interactive(inspect_width,index=widgets.IntSlider(value=25,min=0,max=len(rows)-1,description='Width row')))
else:
    print('Saved reference: optional Jupyter widgets disabled; use the browser library explorer or enable them locally.')

Saved reference: optional Jupyter widgets disabled; use the browser library explorer or enable them locally.

چه چیزی اثبات شد و چه باقی است؟¶

محاسبهٔ تمام‌موج TE، کنترل مستقل لایه، پالایش هر عرض، بازسازی میدان و کتابخانه با مرجع صریح اجرا شده‌اند. scripts/validate_chapter03.py چندمرتبه، آینه، نمونهٔ مرتبهٔ بالاتر و اختیاری grcwa را می‌آزماید. data/chapter-03/validation.json سابقه و محیط آن را دارد.

تمرین‌ها: تفاوت $|t|^2,T$ برای $n_i\ne n_t$ ؛ باقیمانده و خطای مختلط $M=0$ ؛ ارتفاع تازه و پوشش؛ تابش تازه و شمارش و همگرایی؛ ابرسلول همسایهٔ متفاوت با حل‌گر مناسب برای آزمون تقریب موضعی (اینجا نتیجهٔ آن ادعا نمی‌شود)؛ و گروه‌بندی هندسهٔ آموزش و آزمون پیش از فصل ۱۰.

منابع: Moharam و Gaylord، قرارداد grcwa و Meep. حل‌گر دوره اصیل و TE محدود است، نه بستهٔ grcwa. توافق دو RCWA، مقایسهٔ FEM/FDTD یا تجربه نیست.

References¶

Moharam, M. G., & Gaylord, T. K. (1981). Rigorous coupled-wave analysis of planar-grating diffraction. Journal of the Optical Society of America, 71(7), 811–818. 10.1364/JOSA.71.000811