PLaSK library

Loading...
Searching...
No Matches
freecarrier.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include "freecarrier.hpp"
15#include "fd.hpp"
16#include "gauss_matrix.hpp"
17
18 namespace plask { namespace gain { namespace freecarrier {
19
20 constexpr double DIFF_STEP = 0.001;
21
22 OmpNestedLock gain_omp_lock;
23
24 template <typename BaseT>
25 FreeCarrierGainSolver<BaseT>::Level::Level(double E, const Tensor2<double>& M, WhichLevel which, const ActiveRegionParams& params)
26 : E(E), M(M) {
27 thickness = 0.;
28 if (which == EL) {
29 for (size_t i = 0; i < params.U[EL].size(); ++i)
30 if (params.U[EL][i] < E) thickness += params.region.thicknesses[i];
31 } else {
32 for (size_t i = 0; i < params.U[which].size(); ++i)
33 if (params.U[which][i] > E) thickness += params.region.thicknesses[i];
34 }
35}
36
37 template <typename BaseT>
38 FreeCarrierGainSolver<BaseT>::FreeCarrierGainSolver(const std::string& name)
39 : BaseT(name),
40 outGain(this, &FreeCarrierGainSolver<BaseT>::getGainData),
41 outEnergyLevels(this, &FreeCarrierGainSolver<BaseT>::getEnergyLevels),
42 lifetime(0.1),
43 matrixelem(0.),
44 T0(300.),
45 levelsep(0.0001),
46 strained(false),
47 quick_levels(true) {
48 inTemperature = 300.;
49 inTemperature.changedConnectMethod(this, &FreeCarrierGainSolver<BaseT>::onInputChange);
50 inCarriersConcentration.changedConnectMethod(this, &FreeCarrierGainSolver<BaseT>::onInputChange);
51}
52
53 template <typename BaseT> FreeCarrierGainSolver<BaseT>::~FreeCarrierGainSolver() {
54 inTemperature.changedDisconnectMethod(this, &FreeCarrierGainSolver<BaseT>::onInputChange);
55 inCarriersConcentration.changedDisconnectMethod(this, &FreeCarrierGainSolver<BaseT>::onInputChange);
56}
57
58 template <typename BaseT> void FreeCarrierGainSolver<BaseT>::loadConfiguration(XMLReader& reader, Manager& manager) {
59 // Load a configuration parameter from XML.
60 while (reader.requireTagOrEnd()) {
61 std::string param = reader.getNodeName();
62 if (param == "config") {
63 lifetime = reader.getAttribute<double>("lifetime", lifetime);
64 matrixelem = reader.getAttribute<double>("matrix-elem", matrixelem);
65 T0 = reader.getAttribute<double>("T0", T0);
66 strained = reader.getAttribute<bool>("strained", strained);
67 // quick_levels = reader.getAttribute<bool>("quick-levels", quick_levels);
68 if (reader.hasAttribute("substrate")) {
69 substrateMaterial = MaterialsDB::getDefault().get(reader.requireAttribute<std::string>("substrate"));
70 explicitSubstrate = true;
71 }
72 reader.requireTagEnd();
73 } else
74 this->parseStandardConfiguration(reader, manager, "<geometry>, <mesh>, <levels>, or <config>");
75 }
76}
77
78 template <typename BaseT> void FreeCarrierGainSolver<BaseT>::onInitialize() {
79 if (!this->geometry) throw NoGeometryException(this->getId());
80 detectActiveRegions();
81 estimateLevels();
82 outGain.fireChanged();
83}
84
85 template <typename BaseT> void FreeCarrierGainSolver<BaseT>::onInvalidate() {
86 params0.clear();
87 regions.clear();
88 substrateMaterial.reset();
89}
90
91// template <typename BaseT>
92// void FreeCarrierGainSolver<BaseT>::compute()
93//{
94// this->initCalculation(); // This must be called before any calculation!
95// }
96
97 template <typename BaseT>
98 void FreeCarrierGainSolver<BaseT>::ActiveRegionInfo::summarize(const FreeCarrierGainSolver<BaseT>* solver) {
99 holes = BOTH_HOLES;
100 auto bbox = layers->getBoundingBox();
101 total = bbox.upper.vert() - bbox.lower.vert() - bottom - top;
102 materials.clear();
103 materials.reserve(layers->children.size());
104 thicknesses.clear();
105 thicknesses.reserve(layers->children.size());
106 for (const auto& layer : layers->children) {
107 auto block =
108 static_cast<Block<GeometryType::DIM>*>(static_cast<Translation<GeometryType::DIM>*>(layer.get())->getChild().get());
109 auto material = block->singleMaterial();
110 if (!material) throw plask::Exception("{}: Active region can consist only of solid layers", solver->getId());
111 auto bbox = static_cast<GeometryObjectD<GeometryType::DIM>*>(layer.get())->getBoundingBox();
112 double thck = bbox.upper.vert() - bbox.lower.vert();
113 materials.push_back(material);
114 thicknesses.push_back(thck);
115 }
116 double substra = solver->strained ? solver->substrateMaterial->lattC(solver->T0, 'a') : 0.;
117 if (materials.size() > 2) {
118 Material* material = materials[0].get();
119 double e;
120 if (solver->strained) {
121 double latt = material->lattC(solver->T0, 'a');
122 e = (substra - latt) / latt;
123 } else
124 e = 0.;
125 double el0 = material->CB(solver->T0, e, 'G'), hh0 = material->VB(solver->T0, e, 'G', 'H'),
126 lh0 = material->VB(solver->T0, e, 'G', 'L');
127 material = materials[1].get();
128 if (solver->strained) {
129 double latt = material->lattC(solver->T0, 'a');
130 e = (substra - latt) / latt;
131 } else
132 e = 0.;
133 double el1 = material->CB(solver->T0, e, 'G'), hh1 = material->VB(solver->T0, e, 'G', 'H'),
134 lh1 = material->VB(solver->T0, e, 'G', 'L');
135 for (size_t i = 2; i < materials.size(); ++i) {
136 material = materials[i].get();
137 if (solver->strained) {
138 double latt = material->lattC(solver->T0, 'a');
139 e = (substra - latt) / latt;
140 } else
141 e = 0.;
142 double el2 = material->CB(solver->T0, e, 'G');
143 double hh2 = material->VB(solver->T0, e, 'G', 'H');
144 double lh2 = material->VB(solver->T0, e, 'G', 'L');
145 if ((el0 < el1 && el1 > el2) || (hh0 > hh1 && hh1 < hh2) || (lh0 > lh1 && lh1 < lh2)) {
146 if (i != 2 && i != materials.size() - 1) {
147 bool eb = (el0 < el1 && el1 > el2);
148 if (eb != (hh0 > hh1 && hh1 < hh2)) holes = ConsideredHoles(holes & ~HEAVY_HOLES);
149 if (eb != (lh0 > lh1 && lh1 < lh2)) holes = ConsideredHoles(holes & ~LIGHT_HOLES);
150 }
151 if (holes == NO_HOLES)
152 throw Exception("{0}: Quantum wells in conduction band do not coincide with wells is valence band",
153 solver->getId());
154 if ((el0 < el1 && el1 > el2) || (hh0 > hh1 && hh1 < hh2 && holes & HEAVY_HOLES) ||
155 (lh0 > lh1 && lh1 < lh2 && holes & LIGHT_HOLES))
156 wells.push_back(i - 1);
157 } else if (i == 2)
158 wells.push_back(0);
159 if (el2 != el1) {
160 el0 = el1;
161 el1 = el2;
162 }
163 if (hh2 != hh1) {
164 hh0 = hh1;
165 hh1 = hh2;
166 }
167 if (lh2 != lh1) {
168 lh0 = lh1;
169 lh1 = lh2;
170 }
171 }
172 }
173 if (wells.back() < materials.size() - 2) wells.push_back(materials.size() - 1);
174 totalqw = 0.;
175 for (size_t i = 0; i < thicknesses.size(); ++i)
176 if (isQW(i)) totalqw += thicknesses[i];
177
178#ifndef NDEBUG
179 solver->writelog(LOG_DEBUG, "Active region @ {0} ({1}/{2})", origin, bottom, top);
180 assert(materials.size() == thicknesses.size());
181 std::string ws = "-";
182 for (size_t i = 0; i != materials.size(); ++i) {
183 auto w = std::find(wells.begin(), wells.end(), i);
184 if (w != wells.end()) ws = format("{:d}", size_t(w - wells.begin()));
185 solver->writelog(LOG_DEBUG, "[{3}] {0:.2f}nm {1}{2}", thicknesses[i] * 1e3, materials[i]->name(),
186 isQW(i) ? " (QW)" : "", ws);
187 }
188#endif
189}
190
191 template <typename BaseT>
192 FreeCarrierGainSolver<BaseT>::ActiveRegionParams::ActiveRegionParams(const FreeCarrierGainSolver* solver,
193 const ActiveRegionInfo& region,
194 double T,
195 bool quiet,
196 double mt)
197 : region(region) {
198 size_t n = region.materials.size();
199 U[EL].reserve(n);
200 U[HH].reserve(n);
201 U[LH].reserve(n);
202 M[EL].reserve(n);
203 M[HH].reserve(n);
204 M[LH].reserve(n);
205 double substra = solver->strained ? solver->substrateMaterial->lattC(T, 'a') : 0.;
206 Eg = std::numeric_limits<double>::max();
207
208 size_t mi = 0;
209 double me = 0;
210
211 if (!solver->inBandEdges.hasProvider()) {
212 if (!quiet) solver->writelog(LOG_DETAIL, "Band edges taken from material database");
213 size_t i = 0;
214 for (auto material : region.materials) {
215 OmpLockGuard lockq = material->lock();
216 double e;
217 if (solver->strained) {
218 double latt = material->lattC(T, 'a');
219 e = (substra - latt) / latt;
220 } else
221 e = 0.;
222 double uel = material->CB(T, e, 'G'), uhh = material->VB(T, e, 'G', 'H');
223 U[EL].push_back(uel);
224 U[HH].push_back(uhh);
225 double eg = uel - uhh;
226 if (eg < Eg) {
227 Eg = eg;
228 mi = i;
229 me = e;
230 }
231 U[LH].push_back(material->VB(T, e, 'G', 'L'));
232 M[EL].push_back(material->Me(T, e));
233 M[HH].push_back(material->Mhh(T, e));
234 M[LH].push_back(material->Mlh(T, e));
235 ++i;
236 }
237 } else {
238 if (!quiet) solver->writelog(LOG_DETAIL, "Band edges taken from inBandEdges receiver");
239 shared_ptr<Translation<GeometryType::DIM>> shifted(new Translation<GeometryType::DIM>(region.layers, region.origin));
240 auto mesh = makeGeometryGrid(shifted)->getElementMesh();
241 assert(mesh->size() == mesh->vert()->size());
242 auto CB = solver->inBandEdges(BandEdges::CONDUCTION, mesh);
243 auto VB_H = solver->inBandEdges(BandEdges::VALENCE_HEAVY, mesh);
244 auto VB_L = solver->inBandEdges(BandEdges::VALENCE_LIGHT, mesh);
245 for (size_t i = 0; i != mesh->vert()->size(); ++i) {
246 auto material = region.materials[i];
247 OmpLockGuard lockq = material->lock();
248 double e;
249 if (solver->strained) {
250 double latt = material->lattC(T, 'a');
251 e = (substra - latt) / latt;
252 } else
253 e = 0.;
254 double uel = CB[i];
255 double uhh = VB_H[i];
256 U[EL].push_back(uel);
257 U[HH].push_back(uhh);
258 double eg = uel - uhh;
259 if (eg < Eg) {
260 Eg = eg;
261 mi = i;
262 me = e;
263 }
264 U[LH].push_back(VB_L[i]);
265 M[EL].push_back(material->Me(T, e));
266 M[HH].push_back(material->Mhh(T, e));
267 M[LH].push_back(material->Mlh(T, e));
268 }
269 }
270 if (mt == 0.) {
271 if (solver->matrixelem != 0.) {
272 Mt = solver->matrixelem;
273 } else {
274 double deltaSO = region.materials[mi]->Dso(T, me);
275 Mt = (1. / M[EL][mi].c11 - 1.) * (Eg + deltaSO) * Eg / (Eg + 0.666666666666667 * deltaSO) / 2.;
276 if (!quiet) solver->writelog(LOG_DETAIL, "Estimated momentum matrix element to {:.2f} eV m0", Mt);
277 }
278 } else {
279 Mt = mt;
280 }
281}
282
283 template <typename BaseT>
284 double FreeCarrierGainSolver<BaseT>::level(WhichLevel which, double E, const ActiveRegionParams& params, size_t start, size_t stop)
285 const {
286 size_t nA = 2 * (stop - start + 1);
287
288 DgbMatrix A(nA);
289
290 constexpr double fact = 2e-12 * phys::me / (phys::hb_eV * phys::hb_J);
291
292 double m1 = params.M[which][start].c11;
293 double k1_2 = fact * m1 * (E - params.U[which][start]);
294 if (which != EL) k1_2 = -k1_2;
295 double k1 = sqrt(abs(k1_2));
296
297 // Wave functions are confined, so we can assume exponentially decreasing relation in the outer layers
298 A(0, 0) = A(nA - 1, nA - 1) = 1.;
299 A(0, 1) = A(nA - 1, nA - 2) = 0.;
300
301 for (size_t i = start, o = 1; i < stop; i++, o += 2) {
302 double k0_2 = k1_2, k0 = k1, m0 = m1;
303 double d = (o == 1) ? 0. : params.region.thicknesses[i];
304 if (k0_2 >= 0.) {
305 double coskd = cos(k0 * d), sinkd = sin(k0 * d);
306 A(o, o - 1) = coskd;
307 A(o + 1, o - 1) = -sinkd;
308 A(o, o) = sinkd;
309 A(o + 1, o) = coskd;
310 } else {
311 double phi = exp(-k0 * d);
312 A(o, o - 1) = phi;
313 A(o + 1, o - 1) = -phi;
314 A(o, o) = 1. / phi;
315 A(o + 1, o) = 1. / phi;
316 }
317
318 A(o + 2, o) = 0.;
319 A(o - 1, o + 1) = 0.;
320
321 m1 = params.M[which][i + 1].c11;
322 k1_2 = fact * m1 * (E - params.U[which][i + 1]);
323 if (which != EL) k1_2 = -k1_2;
324 if (k1_2 >= 0.) {
325 k1 = sqrt(k1_2);
326 A(o, o + 1) = -1.;
327 A(o + 1, o + 1) = 0.;
328 A(o, o + 2) = 0.;
329 A(o + 1, o + 2) = -(k1 * m0) / (k0 * m1);
330 } else {
331 k1 = sqrt(-k1_2);
332 double f = (k1 * m0) / (k0 * m1);
333 A(o, o + 1) = -1.;
334 A(o + 1, o + 1) = f;
335 A(o, o + 2) = -1.;
336 A(o + 1, o + 2) = -f;
337 }
338 }
339
340 return A.determinant();
341}
342
343 template <typename BaseT>
344 void FreeCarrierGainSolver<BaseT>::estimateWellLevels(WhichLevel which, ActiveRegionParams& params, size_t qw) const {
345 if (params.U[which].size() < 3) return;
346
347 size_t start = params.region.wells[qw], stop = params.region.wells[qw + 1];
348 double umin = std::numeric_limits<double>::max(), umax = std::numeric_limits<double>::lowest();
349 double num = 0.;
350 double ustart, ustop;
351 Tensor2<double> M;
352 for (size_t i = start; i <= stop; ++i) {
353 double ub = params.U[which][i];
354 if (i == start) ustart = ub;
355 if (i == stop) ustop = ub;
356 auto m = params.M[which][i];
357 if (which == EL) {
358 if (ub < umin) {
359 umin = ub;
360 M = m;
361 }
362 } else {
363 if (ub > umax) {
364 umax = ub;
365 M = m;
366 }
367 }
368 if (i != start && i != stop) {
369 double no = 1e-6 / PI * params.region.thicknesses[i] * sqrt(2. * phys::me / (phys::hb_eV * phys::hb_J) * m.c11);
370 num = max(no, num);
371 }
372 }
373 if (which == EL)
374 umax = min(ustart, ustop);
375 else
376 umin = max(ustart, ustop);
377 if (umax < umin)
378 throw Exception("{}: Outer layers of active region have wrong band offset", this->getId()); // TODO make clearer
379 num = 2. * ceil(sqrt(umax - umin) * num); // 2.* is the simplest way to ensure that all levels are found
380 umin += 0.5 * levelsep;
381 umax -= 0.5 * levelsep;
382 double step = (umax - umin) / num;
383 size_t n = size_t(num);
384 double a, b = umin;
385 double fa, fb = level(which, b, params, qw);
386 if (fb == 0.) {
387 params.levels[which].emplace_back(fb, M, which, params);
388 b += levelsep;
389 fb = level(which, b, params, qw);
390 }
391 for (size_t i = 0; i < n; ++i) {
392 a = b;
393 fa = fb;
394 b = a + step;
395 fb = level(which, b, params, qw);
396 if (fb == 0.) {
397 params.levels[which].emplace_back(fb, M, which, params);
398 continue;
399 }
400 if ((fa < 0.) != (fb < 0.)) {
401 boost::uintmax_t iters = 1000;
402 double xa, xb;
403 std::tie(xa, xb) = toms748_solve([&](double x) { return level(which, x, params, qw); }, a, b, fa, fb,
404 [this](double l, double r) { return r - l < levelsep; }, iters);
405 if (xb - xa > levelsep) throw ComputationError(this->getId(), "could not find level estimate in quantum well");
406 params.levels[which].emplace_back(0.5 * (xa + xb), M, which, params);
407 }
408 }
409}
410
411 template <typename BaseT>
412 void FreeCarrierGainSolver<BaseT>::estimateAboveLevels(WhichLevel which, ActiveRegionParams& params) const {
413 if (params.U[which].size() < 5) return; // This makes sense with at least two quantum wells
414
416 size_t N = params.U[EL].size() - 1;
417 double umin = std::numeric_limits<double>::max(), umax = std::numeric_limits<double>::lowest();
418 if (which == EL)
419 umax = min(params.U[EL][0], params.U[EL][N]);
420 else
421 umin = max(params.U[which][0], params.U[which][params.U[which].size() - 1]);
422 Tensor2<double> M;
423 for (size_t i : params.region.wells) {
424 if (i == 0 || i == N) continue;
425 double ub = params.U[which][i];
426 if (which == EL) {
427 if (ub < umin) {
428 umin = ub;
429 M = params.M[which][i];
430 }
431 } else {
432 if (ub > umax) {
433 umax = ub;
434 M = params.M[which][i];
435 }
436 }
437 }
438
439 if (umax <= umin) return;
440
441 double num =
442 2. * ceil(1e-6 / PI * params.region.total * sqrt(2. * (umax - umin) * phys::me / (phys::hb_eV * phys::hb_J) * M.c11));
443 umin += 0.5 * levelsep;
444 umax -= 0.5 * levelsep;
445 double step = (umax - umin) / num;
446 size_t n = size_t(num);
447 double a, b = umin;
448 double fa, fb = level(which, b, params);
449 if (fb == 0.) {
450 params.levels[which].emplace_back(fb, M, which, params);
451 b += levelsep;
452 fb = level(which, b, params);
453 }
454 for (size_t i = 0; i < n; ++i) {
455 a = b;
456 fa = fb;
457 b = a + step;
458 fb = level(which, b, params);
459 if (fb == 0.) {
460 params.levels[which].emplace_back(fb, M, which, params);
461 continue;
462 }
463 if ((fa < 0.) != (fb < 0.)) {
464 boost::uintmax_t iters = 1000;
465 double xa, xb;
466 std::tie(xa, xb) = toms748_solve([&](double x) { return level(which, x, params); }, a, b, fa, fb,
467 [this](double l, double r) { return r - l < levelsep; }, iters);
468 if (xb - xa > levelsep) throw ComputationError(this->getId(), "could not find level estimate above quantum wells");
469 params.levels[which].emplace_back(0.5 * (xa + xb), M, which, params);
470 }
471 }
472}
473
474 template <typename BaseT> void FreeCarrierGainSolver<BaseT>::estimateLevels() {
475 params0.clear();
476 params0.reserve(regions.size());
477
478 size_t reg = 0;
479 for (const ActiveRegionInfo& region : regions) {
480 params0.emplace_back(this, region);
481 ActiveRegionParams& params = params0.back();
482 for (size_t qw = 0; qw < region.wells.size() - 1; ++qw) {
483 estimateWellLevels(EL, params, qw);
484 if (region.holes & ActiveRegionInfo::HEAVY_HOLES)
485 estimateWellLevels(HH, params, qw);
486 else
487 params.levels[HH].clear();
488 if (region.holes & ActiveRegionInfo::LIGHT_HOLES)
489 estimateWellLevels(LH, params, qw);
490 else
491 params.levels[LH].clear();
492 }
493 std::sort(params.levels[EL].begin(), params.levels[EL].end(), [](const Level& a, const Level& b) { return a.E < b.E; });
494 std::sort(params.levels[HH].begin(), params.levels[HH].end(), [](const Level& a, const Level& b) { return a.E > b.E; });
495 std::sort(params.levels[LH].begin(), params.levels[LH].end(), [](const Level& a, const Level& b) { return a.E > b.E; });
496 params.nhh = std::min(params.levels[EL].size(), params.levels[HH].size());
497 params.nlh = std::min(params.levels[EL].size(), params.levels[LH].size());
498 estimateAboveLevels(EL, params);
499 estimateAboveLevels(HH, params);
500 estimateAboveLevels(LH, params);
501
502 if (maxLoglevel > LOG_DETAIL) {
503 {
504 std::stringstream str;
505 std::string sep = "";
506 for (auto l : params.levels[EL]) {
507 str << sep << format("{:.4f}", l.E);
508 sep = ", ";
509 }
510 this->writelog(LOG_DETAIL, "Estimated electron levels for active region {:d} (eV): {}", reg++, str.str());
511 }
512 {
513 std::stringstream str;
514 std::string sep = "";
515 for (auto l : params.levels[HH]) {
516 str << sep << format("{:.4f}", l.E);
517 sep = ", ";
518 }
519 this->writelog(LOG_DETAIL, "Estimated heavy hole levels for active region {:d} (eV): {}", reg - 1, str.str());
520 }
521 {
522 std::stringstream str;
523 std::string sep = "";
524 for (auto l : params.levels[LH]) {
525 str << sep << format("{:.4f}", l.E);
526 sep = ", ";
527 }
528 this->writelog(LOG_DETAIL, "Estimated light hole levels for active region {:d} (eV): {}", reg - 1, str.str());
529 }
530 }
531
532 if (params.levels[EL].empty()) throw Exception("{}: No electron levels found", this->getId());
533 if (params.levels[HH].empty() && params.levels[LH].empty()) throw Exception("{}: No hole levels found", this->getId());
534 }
535}
536
537 template <typename BaseT> double FreeCarrierGainSolver<BaseT>::getN(double F, double T, const ActiveRegionParams& params) const {
538 size_t n = params.levels[EL].size();
539 const double kT = phys::kB_eV * T;
540 constexpr double fact = phys::me * phys::kB_eV / (2. * PI * phys::hb_eV * phys::hb_J); // 1/µm (1e6) -> 1/cm³ (1e-6)
541
542 double N = 2e-6 * pow(fact * T * params.sideM(EL).c00, 1.5) * fermiDiracHalf((F - params.sideU(EL)) / kT);
543
544 for (size_t i = 0; i < n; ++i) {
545 double M = params.levels[EL][i].M.c00;
546 N += 2. * fact * T * M / params.levels[EL][i].thickness * log(1 + exp((F - params.levels[EL][i].E) / kT));
547 }
548
549 return N;
550}
551
552 template <typename BaseT> double FreeCarrierGainSolver<BaseT>::getP(double F, double T, const ActiveRegionParams& params) const {
553 size_t nh = params.levels[HH].size(), nl = params.levels[LH].size();
554 const double kT = phys::kB_eV * T;
555 constexpr double fact = phys::me * phys::kB_eV / (2. * PI * phys::hb_eV * phys::hb_J); // 1/µm (1e6) -> 1/cm³ (1e-6)
556
557 // Get parameters for outer layers
558 double N = 2e-6 * (pow(fact * T * params.sideM(HH).c00, 1.5) * fermiDiracHalf((params.sideU(HH) - F) / kT) +
559 pow(fact * T * params.sideM(LH).c00, 1.5) * fermiDiracHalf((params.sideU(LH) - F) / kT));
560
561 for (size_t i = 0; i < nh; ++i) {
562 double M = params.levels[HH][i].M.c00;
563 N += 2. * fact * T * M / params.levels[HH][i].thickness * log(1 + exp((params.levels[HH][i].E - F) / kT));
564 }
565
566 for (size_t i = 0; i < nl; ++i) {
567 double M = params.levels[LH][i].M.c00;
568 N += 2. * fact * T * M / params.levels[LH][i].thickness * log(1 + exp((params.levels[LH][i].E - F) / kT));
569 }
570
571 return N;
572}
573
574 template <typename BaseT>
575 void FreeCarrierGainSolver<BaseT>::findFermiLevels(double& Fc, double& Fv, double n, double T, const ActiveRegionParams& params)
576 const {
577 double Ue = params.sideU(EL), Uh = params.sideU(HH);
578 double fs = 0.05 * abs(Ue - Uh);
579 if (fs <= levelsep) fs = 2. * levelsep;
580 if (isnan(Fc)) Fc = Ue;
581 if (isnan(Fv)) Fv = Uh;
582 boost::uintmax_t iters;
583 double xa, xb;
584
585 iters = 1000;
586 std::tie(xa, xb) = fermi_bracket_and_solve([this, T, n, &params](double x) { return getN(x, T, params) - n; }, Fc, fs, iters);
587 if (xb - xa > levelsep) throw ComputationError(this->getId(), "could not find quasi-Fermi level for electrons");
588 Fc = 0.5 * (xa + xb);
589
590 iters = 1000;
591 std::tie(xa, xb) = fermi_bracket_and_solve([this, T, n, &params](double x) { return getP(x, T, params) - n; }, Fv, fs, iters);
592 if (xb - xa > levelsep) throw ComputationError(this->getId(), "could not find quasi-Fermi level for holes");
593 Fv = 0.5 * (xa + xb);
594}
595
596 template <typename BaseT>
597 Tensor2<double> FreeCarrierGainSolver<BaseT>::getGain0(double hw,
598 double Fc,
599 double Fv,
600 double T,
601 double nr,
602 const ActiveRegionParams& params) const {
603 constexpr double fac = 1e4 * phys::qe * phys::qe / (2. * phys::c * phys::epsilon0 * phys::hb_J); // 1e4: 1/µm -> 1/cm
604 const double ikT = (1. / phys::kB_eV) / T;
605 const double Dhw = hw - params.Eg;
606
607 Tensor2<double> g(0., 0.);
608
609 for (size_t i = 0; i < params.nhh; ++i) {
610 const double Ec = params.levels[EL][i].E, Ev = params.levels[HH][i].E;
611 const double Ep = hw - (Ec - Ev);
612 if (Ep < 0.) continue;
613 const double sin2 = (Dhw > 0.) ? Ep / Dhw : 0.;
614 const Tensor2<double> pp(1. - 0.5 * sin2, sin2);
615 const double mu = 1. / (1. / params.levels[EL][i].M.c00 + 1. / params.levels[HH][i].M.c00);
616 const double Ecp = Ec + Ep * mu / params.levels[EL][i].M.c00, Evp = Ev - Ep * mu / params.levels[HH][i].M.c00;
617 g += mu * (1. / (exp(ikT * (Ecp - Fc)) + 1) - 1. / (exp(ikT * (Evp - Fv)) + 1)) * pp;
618 }
619
620 for (size_t i = 0; i < params.nlh; ++i) {
621 const double Ec = params.levels[EL][i].E, Ev = params.levels[LH][i].E;
622 const double Ep = hw - (Ec - Ev);
623 if (Ep < 0.) continue;
624 const double sin2 = (Dhw > 0.) ? Ep / Dhw : 0.;
625 const Tensor2<double> pp(0.3333333333333333333333 + 0.5 * sin2, 1.3333333333333333333333 - sin2);
626 const double mu = 1. / (1. / params.levels[EL][i].M.c00 + 1. / params.levels[LH][i].M.c00);
627 const double Ecp = Ec + Ep * mu / params.levels[EL][i].M.c00, Evp = Ev - Ep * mu / params.levels[LH][i].M.c00;
628 g += mu * (1. / (exp(ikT * (Ecp - Fc)) + 1) - 1. / (exp(ikT * (Evp - Fv)) + 1)) * pp;
629 }
630 return fac / (hw * nr * params.region.totalqw) * params.Mt * g;
631}
632
633 template <typename BaseT>
634 Tensor2<double> FreeCarrierGainSolver<BaseT>::getGain(double hw,
635 double Fc,
636 double Fv,
637 double T,
638 double nr,
639 const ActiveRegionParams& params) const {
640 if (lifetime == 0) return getGain0(hw, Fc, Fv, T, nr, params);
641
642 const double E0 = params.levels[EL][0].E - ((params.region.holes == ActiveRegionInfo::BOTH_HOLES)
643 ? std::max(params.levels[HH][0].E, params.levels[LH][0].E)
644 : (params.region.holes == ActiveRegionInfo::HEAVY_HOLES) ? params.levels[HH][0].E
645 : params.levels[LH][0].E);
646
647 const double b = 1e12 * phys::hb_eV / lifetime;
648 const double tmax = 32. * b;
649 const double tmin = std::max(-tmax, E0 - hw);
650 double dt = (tmax - tmin) / 1024.; // TODO Estimate integral precision and maybe chose better integration
651
652 Tensor2<double> g = Tensor2<double>(0., 0.);
653 for (double t = tmin; t <= tmax; t += dt) {
654 // L(t) = b / (π (t²+b²)),
655 g += getGain0(hw + t, Fc, Fv, T, nr, params) / (t * t + b * b);
656 }
657 g *= b * dt / PI;
658
659 return g;
660}
661
662 template struct PLASK_SOLVER_API FreeCarrierGainSolver<SolverWithMesh<Geometry2DCartesian, MeshAxis>>;
663 template struct PLASK_SOLVER_API FreeCarrierGainSolver<SolverWithMesh<Geometry2DCylindrical, MeshAxis>>;
664 template struct PLASK_SOLVER_API FreeCarrierGainSolver<SolverOver<Geometry3D>>;
665
666}}} // namespace plask::gain::freecarrier