Lowest Order Gauss-Mode The unique solutions of ∂ q ∂ z = 1 and ∂ A ∂ z = − A · 1 q are q ( z ) = q 0 + z , where q 0 and z 0 are constants. q 0 q ( z ) . A ( z ) = A 0 Thus, lowest order Gauss mode is e − ikz Ψ ( x, y, z ) E ( x, y, z ) = − z − x 2 + y 2 � � �� q 0 = A 0 q 0 + z exp ik 2( q 0 + z ) Let us normalize the amplitude of this mode by q 0 A 0 = 1 . Then, z + x 2 + y 2 � � �� 1 E ( x, y, z ) = q 0 + z exp − ik 2( q 0 + z ) . – p.25/116
De fi nition of Spot Size De fi nition 1. The spot size is de fi ned by the radius r such that e − 1 = | E ( z, r ) | | E ( z, 0) | . – p.26/116
Spot Size and Beam Waist z + x 2 + y 2 1 � � �� E ( x, y, z ) = q 0 + z exp − ik 2( q 0 + z ) Write 1 1 1 λ q 0 + z = q ( z ) = R ( z ) − i π w ( z ) where � 2 Im ( q 0 ) 2 w ( z ) � ( Re ( q 0 ) + z ) R ( z ) = 1 + ( Re ( q 0 ) + z ) 2 z Im ( q 0 ) + ( Re ( q 0 ) + z ) 2 � � λ − Re ( q 0 ) w ( z ) = Im ( q 0 ) π � � �� z + x 2 + y 2 Phase shift: exp − ik 2 R ( z ) . – p.27/116
Energy of the Beam A short calculation shows: | q 0 + z | 2 = π | Im ( q 0 ) | | w ( z ) | λ and: R 2 | E | 2 d ( xy ) = | A 0 q 0 | 2 π 2 � π λ 2 k | Im ( q 0 ) | 2 Thus, the energy at a slice z = constant is independent of z . . – p.28/116
Types of Resonators There exists several types of resonators . Here, let us study a one way resonator. Other resonators can be transformed to a one way resonator. Let Ω = Ω 2 × [0 , L ] be a res- onator geometry. Let us assume that there are n apertures in the resonator. The start points of these aper- tures are free space free space free space free space z 0 z 1 z 3 z 5 z 7 z 2 z 4 z 6 0 = z 0 ≤ z 1 ≤ z 2 ≤ ... ≤ z n = L. start lense mirror lense mirror x 2 + y 2 � � �� 1 E i ( x, y, z ) = A i q i + ( z − z i ) exp − ik ( z − z i ) + 2( q i + ( z − z i )) where A i := A i q i . . – p.29/116
ABCD Matrices The change of the Gauss-mode is described by ABCD matrices � � A i B i M i = C i D i Then, the beam parameter q i changes as follows q i = A i q i − 1 + B i C i q i − 1 + D i =: M i [ q i − 1 ] . Lemma 1. M i +1 [ M i [ q i − 1 ]] = ( M i +1 M i )[ q i − 1 ] . – p.30/116
Ray Optics and ABCD Matrices An optical ray can be described by the radius r ( z ) and the slope r ′ ( z ) . The change of an optical ray is described by ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎝ r out ⎝ A B ⎝ r in ⎠ = ⎠ ⎠ r ′ C D r ′ out in r out r in Example 1 (Ray-matrix of free space) . ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ L ⎝ r out ⎝ 1 ⎝ r in ⎠ = n 0 ⎠ ⎠ r ′ r ′ 0 1 L out in n 0 . – p.31/116
ABCD matrix of free space Formula 2 (ABCD matrix of free space) . � � � � A B 1 z i − z i − 1 = C D 0 1 and A i = A i − 1 exp( ik ( − ( z i − z i − 1 ))) . – p.32/116
ABCD Matrix of a Lense Formula 3 (ABCD matrix of a lense) . � � � � 1 0 A B 1 = A i = A i − 1 and − 1 1 − 1 1 C D f q i − 1 f d R 2 R 1 r r n 2 , λ 2 n 1 , λ 1 n 1 , λ 1 s 2 s 1 . – p.33/116
ABCD Matrix of a Mirror Formula 4 (ABCD Matrix of a mirror) . � � � � A B 1 0 = − 2 C D 1 R R r s R . – p.34/116
Other ABCD Matrices Formula 5 (ABCD Matrix of a Duct) . Let k = ω √ µ � n ( x ) , where n ( x ) = n 0 − 1 2 n 2 x 2 . Then ( n 0 γ ) − 1 sin( γ z ) � � � � A B cos( γ z ) = , C D − ( n 0 γ ) sin( γ z ) cos( γ z ) where γ 2 = n 2 /n 0 . . – p.35/116
Ray (or Beam) Matrix of the Resonator Using the ABCD matrix M i of each aperture on can calculate the ABCD matrix of the whole resonator by n � � A B � M = M i =: C D i =1 Lemma 2. � � A B = det ( M ) = 1 det C D Proof. Observe that for every aperture the corresponding ABCD matrix M i satis fi es det ( M i ) = 1 . Let r 0 be a start vector. Consider r s = M s r 0 . – p.36/116
Stability Ray of the Resonator Let q a , q b be the eigenvectors of M . Then, r s = c a λ s a q a + c b λ s b q b . Stable Laser: − 1 ≤ | m | ≤ 1 . Then, r s = e i Θ n c a q a + e − i Θ n c b q b , + − i Θ . where λ a,b = e Unstable Laser: | m | ≥ 1 . Then, r s = M s c a q a + M − s c b q b , √ m 2 − 1 . 1 where M = λ a , M = λ b , M = m + . – p.37/116
Exact Solution in a Gaussian “Duct” The refraction index of a Gaussian duct is : k = k 0 (1 − 1 2 n 2 r 2 ) The paraxial approximation and neglecting the small high order term 2 r 2 leads to 4 n 2 1 ∂ Ψ ∂ z − k 0 n 2 r 2 Ψ = 0 △ xy Ψ − 2 ik 0 An exact solution of this equation is: − x 2 + y 2 � � + i λ z Ψ ( x, y, z ) = exp w 2 w 1 1 where w 2 1 2 1 = 2 k 0 √ n 2 and λ = k 0 . . – p.38/116
The Guoy Phase Shift Let us de fi ne the Guoy phase shift ψ ( z ) by: i | q ( z ) | = exp( i ψ ( z )) . q ( z ) This implies tan ψ ( z ) = π w ( z ) 2 R ( z ) λ . Thus, ψ ( z ) = 0 at the waist of the Gaussian beam. Then, one can show 1 q ( z ) = exp( i ( ψ ( z ) − ψ 0 )) q 0 , w 0 w ( z ) where ψ 0 = ψ (0) and q 0 = q (0) . . – p.39/116
Notation in “Lasers and Electro-Optics” In this book the spot size at the waist z = 0 is: � λ z � � 2 � w 2 w 2 D ( z ) = 1 + 0 π w 2 0 By ( ?? ), we get Im ( q 0 ) + ( Re ( q 0 ) + z ) 2 � � � λ � w 2 D ( z ) = w ( z ) = � � � Re ( q 0 )=0 Im ( q 0 ) � Re ( q 0 )= π 0 = λ w 2 π Im ( q 0 ) ⇒ and � 2 Im ( q 0 ) 2 � R ( z ) = ( Re ( q 0 ) + z ) 1 + ( Re ( q 0 ) + z ) 2 . – p.40/116
Hermite-Gaussian Modes � √ � √ w 0 2 x 2 y � � = w H m H n Ψ m,n w w � 1 � �� w 2 + ik − i ( kz − Φ ) − r 2 exp 2 R where � λ z � ( m + n + 1) tan − 1 Φ ( m, n, z ) = π w 2 0 H 0 ( x ) = 1 , H 1 ( x ) = x, H 2 ( x ) = 4 x 2 − 2 , ... H n ( x ) = ( − 1) n e x 2 d n dx n ( e − x 2 ) n = 0 , 1 , ... The set of these functions forms a basis. . – p.41/116
The Laguerre-Gaussian Modes � l � √ 2 r 2 r 2 � � 2 r L l w 2 D cos( l φ ) | Ψ m,n | = E 0 e p w 2 w D D where r, φ are the angle coordinates and L l L l 0 ( x ) = 1 1 ( x ) = l + 1 − x 2 ( x ) = 1 2( l + 1)( l + 2) − ( l + 2) x + 1 2 x 2 L l L n ( x ) = e x d n dx n ( x n e − x ) n = 0 , 1 , ... The set of these functions forms a basis. . – p.42/116
Thermal lensing The refraction index n c ( x ) of a laser crystal changes by a) thermal lensing . b) internal change of the refraction index caused by deformation c) deformation of the end faces of the laser crystal . – p.43/116
Thermal lensing a) The refraction index of a laser crystal changes by temperature Let T 0 be the temperature before heating (refraction index n 0 ). Let T be the temperature caused by the pumping process (refraction index n ). Let η T be the thermal index gradient. (Example: η T = 2 . 2 · 10 − 6 · ◦ C − 1 for Cr 4+ ). Then, n ( x, y, z ) = n 0 + η T ( T ( x, y, z ) − T 0 ) . – p.44/116
Deformation of a Laser Crystal Let B ⊂ R 3 be the original domain of the laser crystal. Let T : B → R 3 be the mapping of the laser deformation such that � x ∈ B � � � T ( x ) + x is the deformed domain of the laser crystal. Heat and deformation of the crystal lead to a refraction index n c ( x ) , x ∈ B such that k c ( x ) = ω √ µ � n c ( x ) . . – p.45/116
Parabolic Fit Assume that B = D × ]0 , L [ , L length of the laser crystal. b) The parabolic fi t of the refraction index is Subdivide ]0 , L [ in N intervals of meshsize h = L N . Let D h be the discretization grid. For every i = 0 , ..., N − 1 : Find n 0 ,i , n 2 ,i such that: � � � n c ( x, y, h ( i + 1 2)) − ( n 0 ,i − 1 2 n 2 ,i ( x 2 + y 2 )) � � � � � l 2 ( D h ) Each of the parameters n 0 ,i , n 2 ,i lead to a matrix ⎡ ⎤ n 0 γ − 1 cos γ i z sin γ i z i A i = ⎣ ⎦ n 0 γ i sin γ i z cos γ i z c) Additionally, perform a parabolic fi t of T ( x, y, 0) and T ( x, y, L ) . . – p.46/116
Beam Propagation Method BPM The paraxial approximation leads to − ∂ 2 Ψ ∂ x 2 − ∂ 2 Ψ ∂ Ψ ∂ z + ( k 2 0 − k 2 ) Ψ = 0 . ∂ y 2 + 2 ik 0 Let us write this equation as follows: ∂ z = ∂ 2 Ψ ∂ x 2 + ∂ 2 Ψ ∂ Ψ ∂ y 2 − ( k 2 0 − k 2 ) Ψ . 2 ik 0 Let Ω = D × ]0 , L [ , then one can apply FE or FD in x, y -direction Crank-Nicolson in z-direction. . – p.47/116
Beam Propagation Method BPM Let Ψ l ( x, y ) be the approximation of Ψ ( x, y, τ l ) , where τ is the time step. Then, Ψ l ( x, y ) is de fi ned by the equations: Ψ l +1 − Ψ l � ∂ 2 Ψ l +1 + ∂ 2 Ψ l +1 1 + ( k 2 0 − k 2 ) Ψ l +1 + 2 ik 0 = ∂ x 2 ∂ y 2 2 τ ∂ 2 Ψ l + ∂ 2 Ψ l � + ( k 2 0 − k 2 ) Ψ l ∂ x 2 ∂ y 2 Ψ initial ( x, y ) Ψ 0 ( x, y ) (initial condition) = Additional boundary conditions are needed . Lenses and mirrors can be discretized by a phase shift. . – p.48/116
Iteration Method of Fox and Li Let Ψ initial be an initial condition at the left mirror. By the BPMethod calculate the beam con fi guration at the right mirror and the re fl ected beam con fi guration Ψ end := B ( Ψ initial ) at the left mirror. If Ψ initial = Ψ end , then Ψ initial is an eigenvector Ψ eigen of the BPM operator B . The iteration method of Fox and Li is a power iteration method for the eigenvalue problem of the BPM operator B . This means: Ψ 1 = Ψ initial , Ψ s +1 = B ( Ψ initial,s ) Ψ eigen s →∞ Ψ s = lim . – p.49/116
Weak Formulation of the Helmholtz Equati Let V := { v ∈ H 1 ( Ω ) | Γ M = 0 } . Then, −△ u − k 2 u = f u | Γ M = 0 , u · ik + ∂ u ∂ n | Γ R = 0 transforms to: � Problem 1. Find u ∈ V = { v ∈ H 1 ( Ω ) � v | Γ M = 0 } such that � � � ∂ u � v − k 2 u ¯ ∇ u ∇ ¯ v dµ − ik ∂ n ¯ v dµ = f ¯ v dµ for every v ∈ V . Ω Γ Ω . – p.50/116
Weak Formulation of the Helmholtz Equati De fi ne the bilinear form ∂ u � � v − k 2 u ¯ a ( u, v ) = ∇ u ∇ ¯ v dµ − ∂ n ¯ v dµ Ω Γ Then, the week form of the Helmholtz equation is transforms to: � Problem 2. Find u ∈ V = { v ∈ H 1 ( Ω ) � v | Γ M = 0 } such that � � a ( u, v ) = f ¯ v dµ for every v ∈ V . Ω . – p.51/116
Properties of a ( u, v ) : a) The local part of a ( u, v ) is the bilinear form � a loc ( u, v ) = v − k 2 u ¯ ∇ u ∇ ¯ v dµ Ω Let k be constant. Then, a loc is not positive de fi nite, since ⎧ if k 1 > k > 0 ⎪ ⎨ + + a loc ( e − ik 1 z , e − ik 1 z ) = if k 1 = k = 0 ⎪ if k 1 < k < 0 ⎩ . – p.52/116
Properties of a ( u, v ) : − ikz are + b) Let k be constant. Then, the functions e contained in the local kernel of a . This means + − ikz , v ) = 0 for every v ∈ H 1 a ( e 0 ( Ω ) . . – p.53/116
Properties of a ( u, v ) : c) The bilinear form a ( u, v ) is H 1 -coercive. This means that there exist c, C > 0 such that Re ( a ( u, u )) + C ∥ u ∥ 2 L 2 ≥ c ∥ u ∥ 2 ∀ u ∈ H 1 ( Ω ) H 1 . – p.54/116
Properties of a ( u, v ) : d) The problem Find u ∈ V such that for every v ∈ V a ( u, v ) = 0 has the unique solution u = 0 , if k > 0 . . – p.55/116
Properties of a ( u, v ) : d) The problem Find u ∈ V such that for every v ∈ V a ( u, v ) = 0 has the unique solution u = 0 , if k > 0 . . – p.56/116
Boundary Conditions Let Ω ⊂ R d , d = 1 , 2 , 3 be an open d -dimensional open bounded domain. Consider −△ u − k 2 u = 0 The rays exp( ik ⃗ m · x ) are solutions of this equation, where m = 1 . ⃗ . – p.57/116
Boundary Conditions in 1D First, let us consider the 1D case d = 1 and Ω =]0 , 1[ . Then and exp( ikz ) exp( − ikz ) are solutions of − ∂ 2 u ∂ z 2 − k 2 u = 0 . Let us assume that the re fl ection of the ray exp( − ikz ) at the point 0 is α exp( ikz ) . This means we need a boundary condition at 0 with solution u ( x ) = exp( − ikz ) + α exp( ikz ) . A suitable boundary condition is u | z =0 (1 − α ) ik + (1 + α ) ∂ u ∂ z | z =0 = 0 . . – p.58/116
Simple Boundary Conditions Re fl ecting boundary condition: u | z =0 = 0 Non-re fl ecting boundary condition: u | z =0 ik + ∂ u n | z =0 = 0 . ∂⃗ . – p.59/116
Non-Re fl ecting Boundary Condition in 2D,3 Observe that x →−∞ exp( − i ( k + i α ) ⃗ lim m · x ) = 0 , where α > 0 . This leads to the concept: Extend the PDE outside of the domain. Add an adsorbtion coef fi cient α outside of the domain. Put homogenous Dirichlet boundary conditions at a certain distance for away from the boundary. . – p.60/116
Dif fi culties of a Pure FE Discretization One dif fi culty is the large number of discretization grid points which are needed in case of long resonators. Dif fi culties occur, if 1 cm = L >> 5 λ = 10 µm . Then, more than 20 ∗ 1000 = 20000 grid points are needed only in z-direction. The second dif fi culty is that a is not symmetric positive de fi nite and the resulting linear equation system cannot ef fi ciently be solved by multigrid or any other standard iterative solver. There exist several eigenvectors with eigenvalues close to each other. A very accurate discretization of the non-re fl ecting boundary condition is needed. . – p.61/116
Modeling the Wave in a Resonator Let us model the wave E ( x, y, z ) in a one way resonator by the following equations: ∂ u − ∆ u + 2 ik f ∂ z + k s (2 k f − k s ) u = ξ u � � E ( x, y, z ) = exp − i ( k f − ε ) z u ( x, y, z ) 2 ε k f = ξ . – p.62/116
Modeling the Wave in a Resonator Let us assume that Φ ⊂ R 2 is a bounded and connected domain with a piecewise smooth boundary and let Ω = Φ × ]0 , L [ , where L > 0 . Let us subdivide the boundaries of Ω by Γ L := Φ × { L } and Γ r := ∂ Ω \ ( Γ 0 ∪ Γ L ) . Γ 0 := Φ × { 0 } , For reasons of simplicity, let us additionally assume that we choose k f such that exp [ jLk f ] = 1 . (7) . – p.63/116
Two-Wave Ansatz for Resonator Modeling Let us model the resonator by a forward wave E r and the backward wave E l such that E = E r + E l , where each of these waves satisfy the Helmholtz equation . This leads to the eigenvalue problem: ∂ u r ∂ z + ( k 2 f − k 2 ) u r − ∆ u r + 2 jk f = ξ u r , (8) ∂ u l ∂ z + ( k 2 f − k 2 ) u l − ∆ u l − 2 jk f = ξ u l , where E r ( x, y, z ) = exp [ − jk f z ] u r ( x, y, z ) , E l ( x, y, z ) = exp [ − jk f ( L − z )] u l ( x, y, z ) , . – p.64/116
Boundary Conditions for Two-Wave Ansat To satisfy the boundary conditions ( ?? ) and ( ?? ), we need the boundary conditions � u r + u l = 0 , (9) � � Γ 0 ∪ Γ L ∂ u r � n − jC b u r = 0 , (10) � ∂⃗ � Γ r ∂ u l � n − jC b u l = 0 . (11) � ∂⃗ � Γ r � Observe that ( ?? ) is needed to obtain E r + E l = 0 from � � Γ 0 ∪ Γ L � u r + u l = 0 . � � Γ 0 ∪ Γ L To obtain a system of equations with enough equations, we additionally need the boundary condition ∂ u r ∂ z − ∂ u l � = 0 . (12) � ∂ z � Γ 0 ∪ Γ L . – p.65/116
Weak Formulation Let us de fi ne � H 1 = � � ⃗ ( u r , u l ) ∈ H 1 ( Ω ) × H 1 ( Ω ) � � u r + u l Γ 0 = 0 , u r + u l | Γ L = 0 . � � a (( u r , u l ) , ( v r , v l )) = ⃗ � � � ∂ u r � v r + ( k 2 f − k 2 ) u r ¯ = ∇ u r ∇ ¯ v r + 2 jk f ∂ z ¯ v r − jC b u r ¯ v r Ω Γ r � � � ∂ u l � v l + ( k 2 f − k 2 ) u l ¯ + ∇ u l ∇ ¯ v l − 2 jk f ∂ z ¯ v l − jC b u l ¯ v l , Ω Γ r where we assume that k ∈ L ∞ ( Ω ) . Now, the weak formulation is: H 1 and ξ ∈ C such that u = ( u r , u l ) ∈ ⃗ Find ⃗ � v = ( v r , v l ) ∈ ⃗ H 1 . ⃗ a ( ⃗ u, ⃗ v ) = ξ u r v r + u l v l ∀ ⃗ Ω . – p.66/116
Properties of a Lemma 3. Let C b = 0 . Then, ⃗ a ( ⃗ u, ⃗ v ) is symmetric. But ⃗ v ) is not positive de fi nite. a ( ⃗ u, ⃗ . – p.67/116
Trilinear Finite Elements Ω h := { ( ih x , jh y , kh z ) | i, j = − N x , ..., N x and k = 0 , ..., N z } , where we set h = ( h x , h y , h z ) . Furthermore, we obtain the following set of cells τ := { [ ih x , ( i + 1) h x ] × [ ih y , ( i + 1) h y ] × [ ih z , ( i + 1) h z ] | and i, j = − N x , ..., N x − 1 k = − N z , ..., N z − 1 } . Let us de fi ne the space of trilinear fi nite elements by � � V h := u ∈ C ( Ω ) ∀ T ∈ τ : ∃ c 1 , c 2 , c 3 , c 4 , c 5 , c 6 , c 7 , c 8 ∈ C : � � u ( x, y, z ) | T = c 1 + c 2 x + c 3 y + c 4 z + � c 5 xy + c 6 yz + c 7 xy + c 8 xyz . – p.68/116
Two Wave Finite Element Space Let us de fi ne the fi nite element space � � � ⃗ ⊂ ⃗ H 1 � V h := ( u h,r , u h,l ) ∈ V h × V h � u h,r + u h,l Γ 0 = 0 , u h,r + u h,l | Γ L = 0 � � An unstable FE-discretization is: u h ∈ ⃗ Find ⃗ V h such that � ⃗ v h ∈ ⃗ a ( ⃗ u h , ⃗ v h ) = f ⃗ vd ∀ ⃗ V h ⃗ Ω The resulting linear equation system is dif fi cult to solve. . – p.69/116
FE-Theory for Sym. Positive De fi nite S.F. Theorem 2. Let a be a continuous symmetric positive de fi nite sesquilinear form on a Hilbert space V , V h a closed subspace and f ∈ V ′ . Furthermore, let u ∈ V and u h ∈ V h such that a ( u, v ) = f ( v ) ∀ v ∈ V a ( u h , v h ) = f ( v h ) ∀ v h ∈ V h Then, ∥ u − u h ∥ E ≤ inf v h ∈ V h ∥ u − v h ∥ E , where ∥ . ∥ E is the norm corresponding to a . . – p.70/116
FE-Theory for Positive De fi nite S.F. Theorem 3. Let a be a continuous positive de fi nite sesquilinear form on a Hilbert space V , V h a closed subspace of V and f ∈ V ′ . Furthermore, let us assume that a is V -elliptic. This means that there is a constant α > 0 such that | a ( u, u ) | ≥ α ∥ u ∥ 2 ∀ u ∈ V. The continuity of a implies that there is a constant C such that a ( u, v ) ≤ C ∥ u ∥∥ v ∥ ∀ u, v ∈ V. Furthermore, let u ∈ V and u h ∈ V h such that a ( u, v ) = f ( v ) ∀ v ∈ V, a ( u h , v h ) = f ( v h ) ∀ v h ∈ V h . Then, ∥ u − u h ∥ ≤ C v h ∈ V h ∥ u − v h ∥ . inf α . – p.71/116
Streamline-Diffusion Discretization ∂ u − ∆ u + 2 ik f ∂ z + k s (2 k f − k s ) u = f on B . (13) Let us extend the subdivision τ of Ω to a subdivision τ B of B by using the same meshsize. Furthermore, let V h, B be the corresponding fi nite element space of trilinear functions. Discretization: Find u h ∈ V h, B such that � − C b u h ¯ v h + ∂ B � ∂ u h v h + h ρ ∂ + ∇ u h ∇ ¯ v h + 2 ik f ∂ z (¯ ∂ z ¯ v h ) + B v h + h ρ ∂ � v h + h ρ ∂ k s (2 k f − k s ) u h (¯ ∂ z ¯ v h ) d = f (¯ ∂ z ¯ v h ) d T ∀ v h ∈ V h, B . . – p.72/116
Streamline-Diffusion Discretization An stable FE-discretization is: u h ∈ ⃗ Find ⃗ V h such that � ∂ u h , r ∂ ¯ v h , r a ( ⃗ u h , ⃗ v h ) + d ⃗ h ρ 2ik f ∂ z ∂ z Ω � ∂ u h , l ∂ ¯ v h , l + d h ρ 2ik f ∂ z ∂ z Ω � ⃗ = f ⃗ v h d Ω �� � � ∂ ∂ v h ∈ ⃗ + h ρ f r ∂ z ¯ v h,r d − f l ∂ z ¯ v h,l d ∀ ⃗ V h Ω Ω We call this discretization streamline-diffusion discretization. However, there are no streamlines. In case of a convection-diffusion equation, this discretization converges with O ( h 2 ) . . – p.73/116
Ellipticity The sesquilinear form of the streamline-diffusion discretization is a h ( ⃗ u h , ⃗ v h ) = ⃗ a ( ⃗ u h , ⃗ v h ) ⃗ � ∂ u h , r ∂ ¯ � ∂ u h , l ∂ ¯ v h , r v h , l + d + h ρ h ρ 2ik f 2ik f d ∂ z ∂ z ∂ z ∂ z Ω Ω v h ∈ ⃗ V h the following inequality holds: Lemma 4. For every ⃗ 2 � � v h ∂⃗ � � | ⃗ a h ( ⃗ v h , ⃗ v h ) | ≥ hk f ρ . � � ∂ z � � . – p.74/116
A Smoothness Result a h ( ⃗ u h , ⃗ v h ) = ⃗ a ( ⃗ u h , ⃗ v h ) ⃗ � ∂ u h , r ∂ ¯ v h , r � ∂ u h , l ∂ ¯ v h , l + d + h ρ h ρ 2ik f 2ik f d ∂ z ∂ z ∂ z ∂ z Ω Ω H 1 such that: h ∈ ⃗ u c Lemma 5. Let ⃗ � ⃗ v ∈ ⃗ H 1 . u c a h ( ⃗ h , ⃗ v ) = f ⃗ v d ∀ ⃗ ⃗ Ω Then, ∂ 2 ⃗ u c � � C ∥ ⃗ h � � ≤ f ∥ L 2 . � � ∂ z 2 h ρ k f � � 2 . – p.75/116
FE-Theory for Positive De fi nite S.F. Since ⃗ a h satis fi es the Garding inequality, one can prove the following convergence theorem: H 1 such that Theorem 4. Assume ⃗ u = ( u r , u l ) ∈ ⃗ f = ( f r , f l ) ∈ L 2 ( Ω ) 2 . Let ⃗ � v = ( v r , v l ) ∈ ⃗ H 1 . a ( ⃗ u, ⃗ v ) = f r v r + f l v l ∀ ⃗ ⃗ Ω u h = ( u r,h , u l,h ) ∈ ⃗ V h such that and ⃗ � v h = ( v r , v l ) ∈ ⃗ ⃗ a h ( ⃗ u h , ⃗ v h ) = f r v r + f l v l ∀ ⃗ V h . Ω u h converges to ⃗ u . Then, ⃗ . – p.76/116
A Symmetry Consideration Instead of ⃗ a h ( ⃗ u h , ⃗ v h ) := ⃗ a ( ⃗ u h , ⃗ v h ) � ∂ u h , r ∂ ¯ � ∂ u h , l ∂ ¯ v h , r v h , l + h ρ d + h ρ 2ik f 2ik f d ∂ z ∂ z ∂ z ∂ z Ω Ω one can de fi ne a h ( ⃗ u h , ⃗ v h ) := a ( ⃗ u h , ⃗ v h ) ⃗ ⃗ � � ∂ u h , r ∂ ¯ v h , r ∂ u h , l ∂ ¯ v h , l − d − h ρ h ρ 2ik f 2ik f d ∂ z ∂ z ∂ z ∂ z Ω Ω Then, the meaning of u h,r and u h,l changes and the meaning of t and − t in the ansatz E ( x, y, z, t ) = exp( i ω t ) E ( x, y, z ) . . – p.77/116
Streamline-Diffusion Discretization in 1D In 1D the sesquilinear form � 1 ∂ u h v h d + h ρ∂ u h ∂ ¯ v h 2 ik f ∂ z ¯ ∂ z d ∂ z 0 leads to the stencil 1 1 ik f 2 ( − 1 0 1) + ik f h ρ h ( − 1 2 − 1) = ik f ( − 1 1 0) for ρ = 1 2 . This is the FD upwind discretization. An exact solver for the resulting equation system is a Gauss-Seidel relaxation from left to right. . – p.78/116
Iterative Solver Hackbusch’s rule: Consider a singular perturbed problem with parameter � → � 0 . Then, construct an iterative solver such that this solver is an exact solver for � 0 (usually � 0 = 0 ). The transformed one way resonator equation is: − � ∆ u + 2 i ∂ u ∂ z + � k s (2 k f − k s ) u = �ξ u 1 where � = k f . Then, in 1D, the streamline diffusion discretization stencil for � → 0 is i ( − 1 1 0) An exact solver for the corresponding equation system with suitable bound- ary conditions is a relaxation from left to right. Thus, we used a relaxation from left to right as a preconditioner for GMRES. . – p.79/116
Numerical Results Figure 1: Gauss-Mode by FE . – p.80/116
Numerical Results Figure 2: Gauss-Mode by FE . – p.81/116
Modeling Optical Apparatuses Here: One-way resonator with a lense or an interface at the point l 0 with 0 < l 0 < L . Let us write and Ω a = Ψ × [0 , l 0 ] ⊂ B Ω b = Ψ × [ l 0 , L ] ⊂ B . Ψ I = Ψ × { l 0 } . Then, the ansatz E ( x, y, z ) = exp [ − ik f z ] u ( x, y, z ) is not appropriate. Instead, we use the ansatz � � � for z < l 0 exp − ik f, Ω a z u ( x, y, z ) E ( x, y, z ) = u ( x, y, z ) for z > l 0 . , � � exp − ik f, Ω b z u ( x, y, z ) where k f,a is an average value of k f in Ω a and k f,b is an average value of k f in Ω b . . – p.82/116
Modeling Optical Apparatuses � ∂ u a Ξ ( u, v ) := ∇ u ∇ ¯ v + 2 ik f, Ξ ∂ z ¯ v + k s (2 k f, Ξ − k s ) u ¯ v d Ξ a ( u, v ) = a Ω a ( u, v ) + a Ω b ( u, v ) . Phase shift of the apparatuses: ϕ ( x, y ) . Then, let us de fi ne the space � � u ∈ L 2 ( B ) u | Ω a ∈ H 1 ( Ω a ) , u | Ω b ∈ H 1 ( Ω b ) H ab = and � � � u | Ψ I , Ω a · ϕ = u | Ψ I , Ω b . Find u ∈ H ab such that � a ( u, v ) = u ¯ v d ∀ v ∈ H ab . B . – p.83/116
Gain and Absorption To simulate gain and absorption in the Helmholtz equation −△ u − k 2 u = 0 we apply the formula k 2 = ω 2 µ � + j 2 ω √ µ �α . By rate equations, we obtain K = σ c and α gain = 1 / 2 σ N Thus, we get ω 2 µ � + j ω √ µ � ( σ N − 2 α absorption ) k 2 = ω 2 µ � + j ω √ µ � ( σ N − 1 = ) . τ c . – p.84/116
Time-Dependent Behavior Using the ansatz E ( x, y, z, t ) = exp( i ω t )( E r ( x, y, z, t ) + E l ( x, y, z, t )) we obtain µ �∂ 2 u r ∂ t 2 + iµ �ω ∂ u r ∂ u r ∂ z − ( k 2 f − k 2 ) u r , = ∆ u r − 2 jk f ∂ t µ �∂ 2 u l ∂ t 2 + iµ �ω ∂ u l ∂ u l ∂ z − ( k 2 f − k 2 ) u l , = ∆ u l + 2 jk f ∂ t ∂ N − γ Nn σ c − N + N tot ( γ − 1) = + R p ( N tot − N ) ∂ t τ f ω 2 µ � + j ω √ µ � ( σ N − 1 k 2 = ) τ c � 2 � ω | E | 2 n = | u r | 2 + | u l | 2 . | E | 2 = . – p.85/116
Weak Formulation for the Maxwell Equatio The time-periodic vector Helmholtz equation is E − k 2 ⃗ E = ⃗ ∇ × ∇ × ⃗ f. The bilinear form of the weak formulation is positive de fi nite: � E · ∇ × ¯ E ¯ W + k 2 ⃗ a ( ⃗ E, ⃗ ∇ × ⃗ ⃗ ⃗ W ) = W d ( x, y, z ) Ω Then, we obtain � u, e − ikz ⃗ v + ( k 2 − k 2 u · ∇ × ¯ a ( e − ikz ⃗ w ) = ∇ × ⃗ f ) u z ¯ v z ⃗ Ω � ∂ u z ∂ y − ∂ u y � v y + ( k 2 − k 2 − ik f 2 ¯ f ) u y ¯ v y ∂ z � ∂ u z ∂ x − ∂ u x � v x + ( k 2 − k 2 − ik f 2 ¯ f ) u x ¯ v x d ( x, y, z ) ∂ z . – p.86/116
VCSEL (Vertical Cavity Surface Emitting Laser) light output p − contact oxide aperture top DBR passivation active layer bottom DBR current flow substrate n − contact . – p.87/116 M FO OP
DFB Laser (Distributed Feedback Laser) injection current contact L w λ 4 active layer substrate . – p.88/116 Fi 4 VCSEL (Di t ib t d F db k L )
Distributed Bragg Re fl ectors (DBR) Let us assume that the resonator has the form Ω = Ψ × [0 , L ] and that 0 = l 0 < l 1 < ... < l s = L Furthermore, let us assume that the resonator has the refraction index n i ( k i ) in the layer Ψ × [ l i − 1 , l i ] . Assume that − E ′′ − k 2 E = 0 . Let us assume the k is constant in the interior of [ l i − 1 , l i ] . Then, for z ∈ [ l i − 1 , l E ( z ) = c i,r exp( − ik i ( z − l i − 1 ))+ c i,l exp( ik i ( z − l i − 1 )) . – p.89/116
Transmission Matrix of One Layer c i,r c i +1 ,r n i +1 n i c i,l c i +1 ,l ⎛ ⎞ ⎛ ⎞ ⎝ c i +1 ,r ⎝ c i,r ⎠ = M i ⎠ c i +1 ,l c i,l h i ⎛ ⎞ ⎛ ⎞ ⎝ k i +1 + k i k i +1 − k i ⎝ exp( − ik i h i ) 0 1 ⎠ · ⎠ . M i = 2 k i +1 k i +1 − k i k i +1 + k i 0 exp( ik i h i ) . – p.90/116
Transmission Matrix and Scattering Matri c 1 ,r c 2 ,r black box c 1 ,l c 2 ,l In general one can describe the behavior by a scattering matrix S and a transmission matrix T : � � � � � � � � c 1 ,r c 2 ,r c 2 ,r c 1 ,r = T = S c 1 ,l c 2 ,l c 1 ,l c 2 ,l . – p.91/116
Re fl ection Property of DFB Example 2. Let us study 101 layers with refraction index n 0 , n 1 , n 0 , ..., n 0 , √ � 0 µ 0 n 0 , where √ � 0 µ 0 = 1 λ 0 = 1 . 6 · 10 − 6 , k 0 = 2 π k λ 0 , and ω = c and n 0 = 3 . 277 . Let us choose c 2 ,l = 1 , c 1 ,r = 0 . Then, c 1 ,l shows the behavior of the construction. A high re fl ectivity is obtained for ω = ω 0 , 3 ω 0 , 5 ω 0 , ... . Re fl ection behavior for n 1 = Re fl ection behavior for n 1 = 3 . 275 . 3 . 220 . . – p.92/116
Finite Elements of DFB Let Ω h be a grid of meshsize h for the domain Ω = [0 , L ] . Furthermore, let v p be the nodal basis function with respect to linear elements. Then, de fi ne v l e ikz v p = p v r e − ikz v p = p ⎧ e ikz v p ( z ) for z ≤ p and ⎨ v m = p e − ikz v p ( z ) for z > p . ⎩ Now, let us de fi ne the FE space V ref = span { v l p , v r p , v m p | p ∈ Ω h } . h This FE space leads to the results as the transfer matrix method. But these basis functions can be extended to 2D and 3D. . – p.93/116
Time Discretization Let us recall the scalar Helmholtz equation ( ?? ): E = − µ � ∂ 2 � � −△ ˜ ˜ E . ∂ t 2 The ansatz ˜ E ( x, y, z, t ) = exp( i ω t ) E ( x, y, z, t ) leads to µ �∂ 2 E ∂ t 2 + iµ �ω∂ E ∂ t = △ E + µ �ω 2 E. Since ω 2 is large in comparison to µ � , we apply the following model: iµ �ω∂ E ∂ t = △ E + k 2 E. . – p.94/116
Time Discretization Crank-Nicolson discretization of this equation leads to iµ �ω E s +1 − E s = 1 △ E s + k 2 E s + △ E s +1 + k 2 E s +1 � � . 2 τ Let us analyze this equation by Fourier analysis in 2D. Then, for E s = a s sin( l x x ) sin( l y y ) , we obtain y + k 2 ) + i µ �ω 1 2 ( l 2 x + l 2 a s +1 = τ a s y + k 2 ) − i µ �ω 1 2 ( l 2 x + l 2 τ This equation implies | a s +1 | = | a s | if k ∈ R . This means a real k does not lead to a gain or an absorption. An explicit or implicit Euler discretization does not have this property. . – p.95/116
Modeling the Wave in a Resonator Let us assume that Ω = Ω 2 D × [0 , L ] is the domain of a laser resonator, where L is the length of the resonator. Here, let us assume that E 1 , ..., E M are eigenmodes obtained by a Gauss mode analysis or another method. Thus, E i : Ω → C are functions, which we normalize as follows � | E i | 2 d ( x, y, z ) = 1 . Ω . – p.96/116
Model Assumption 1 The electrical fi eld E of the total optical wave is approximated by M eigenmodes: M � E ( t, x, y, z ) = ξ i ( t ) E i ( x, y, z ) , i =1 where ξ i : [ t 0 , ∞ [ → R is the time-dependent coef fi cient of the i -th mode. Then, the photon density of the mode ξ i ( t ) E i ( x, y, z ) is � � | ξ i ( t ) E i ( x, y, z ) | 2 = Ξ i ( t ) | E i ( x, y, z ) | 2 , n i ( t, x, y, z ) = 2 � ω i 2 � ω i where we abbreviate Ξ i ( t ) = | ξ i ( t ) | 2 . . – p.97/116 i th f f th th d
Model Assumption 2 The modes are incoherent modes. Here, this means that the total photon density n ( t, x, y, z ) can be written as M � n ( t, x, y, z ) = n i ( t, x, y, z ) . (14) i =1 . – p.98/116
Model Assumption 3 The local photon densities n i ( t, x, y, z ) and the population inversion density N ( t, x, y, z ) satisfy the rate equations: ∂ n i Nn i σ c − n i = + S, i = 1 , .., M, (15) ∂ t τ c ∂ N − γ Nn σ c − N + N tot ( γ − 1) = + R pump ( N tot − N ) . (16) ∂ t τ f . – p.99/116
ODE System +2 � ω i ∂ Ξ i � N | E i | 2 d ( x, y, z ) σ c − Ξ i � ∂ t = Ξ i S d ( x, y, z ) , i = τ c � Ω Ω (17) M ∂ N Ξ i | E i | 2 − N + N tot ( γ − 1) � � = − γ N σ c + R pump ( N tot − N ∂ t 2 � ω i τ f i =1 (18) These equations form a solvable system of ordinary differ- ential equations, which describes the time-dependent be- havior of M modes. This behavior is mainly in fl uenced by the pump con fi guration R pump . . – p.100/116
Recommend
More recommend
