The selection rules are important for classifying beta decay's transitions. They are based on the conservation of angular momentum and parity during decay.

**The selection rules are important for classifying beta decay's transitions. They are based on the conservation of angular momentum and parity during decay.**

## Beta + vs Beta - Decay

There exist different kinds of beta decay, two of them being beta + ($\beta^+$) and beta - ($\beta^-$) decay. In $\beta^+$ decay, a proton $p^+$ decays into a neutron $n$, a positron $e^+$, and an electron neutrino ${\nu}_e$:

$$
    p^+ \rightarrow n + e^+ + \nu_e.
$$

In $\beta^-$ decay, a neutron $n$ decays into a proton $p^+$, an electron $e^-$, and an electron antineutrino $\bar{\nu}_e$:

$$
    n \rightarrow p + e^- + \bar{\nu}_e.
$$

## Fermi and Gamow-Teller Transitions

For $\beta$ decay, the total spin quantum number $S$ is obtained by summing the spins of the emitted particles. Since both the electron and the neutrino have a spin of 1/2, the sum will therefore be either 0 (antiparallel) or 1 (parallel), depending on the orientation of the spins. $S=0$ is called a Fermi transition, while $S=1$ is called a Gamow-Teller transition. 

If we take a look at the change in angular momentum of the nucleus

$$
    \Delta I = |I_i - I_f|,
$$

where $I_i$ is the angular momentum of the initial nuclear state and $I_f$ is the angular momentum of the final nuclear state. We see that $\Delta I$ does not change ($\Delta I = 0$) for a Fermi transition since $S=0$. When the electron and neutrino spins are parallel (Gamow-Teller decay) they carry a total angular momentum of 1 unit so the conversation of angular looks like this in vector form:

$$
    \vec{I_i} = \vec{I_f} + \vec{1}.
$$

This is only possible if $\Delta I = 0$ or $1$. However, it is not possible when $I_i = I_f = 0$ ($0 \rightarrow 0$), in that case only the Fermi transition can contribute.

## Allowed and Forbidden Transitions

The orbital angular momentum $L_\beta$ can take on all positive integers and also dictates the type of transition. When $L_\beta=0$ the decay is called allowed, while the remaining values are called forbidden (because they are generally less probable than allowed decays), with $L_\beta=1$ being the first-forbidden decay, $L_\beta=2$ being the second-forbidden, and so on. The allowed and all the forbidden decays have each a Fermi and Gamow-Teller type, and the allowed Fermi decay is sometimes called superallowed. 

## Parity in Beta Decay

$L_\beta$ also determines the parity change of the wave function for beta decay according to $\Delta \pi=(-1)^{L_\beta}$. So for example, allowed transitions' parity does not change.

## First-forbidden Decay

First-forbidden decays, $L_\beta = 1$, like allowed decays, have both a Fermi and Gamow-Teller type. For the Fermi transition type, we couple $S=0$ with $L_\beta = 1$ to get a total of 1 unit of angular momentum carried by the beta decay so the change in angular momentum becomes $\Delta I = 0$ or $1$ (but no $0 \rightarrow 0$ transition). For the Gamow-Teller type $S=1$ combined with $L_\beta = 1$ gives a total of 0, 1, or 2 units of total angular momentum so $\Delta I = 0$, $1$, or $2$. Since $L_\beta = 1$ is an odd number the parity changes so $\Delta \pi=(-1)^{L_\beta} = $ yes.

## Second-forbidden Decay

Using a similar methodology as the first-forbidden decay if we couple $L_\beta = 2$ and $S=0$ we in principle get $\Delta I = 0$, $1$, or $2$ for the Fermi type, and for the Gamow-Teller type ($L_\beta = 2$ and $S=1$) we in principle get $\Delta I = 0$, $1$, $2$, or $3$. However, since the parity does not change ($\Delta \pi=(-1)^{L_\beta} = $ no) the $\Delta I = 0$ and $1$ cases fall within the allowed decay's selection rules so we can exclude them since the allowed variations are much more probable to occur. 

## Beta Decay Transition Rules

A table consolidating all the transition rules for beta decay mentioned above is shown below.

![Beta Decay Transitions, Allowed, Forbidden, Fermi, Gamow-Teller, Parity Change](https://i.imgur.com/CPhtzE3.png "Beta Decay Transitions, Allowed, Forbidden, Fermi, Gamow-Teller, Parity Change")

Allowed, Forbidden, Fermi, and Gamow-Teller Beta Decay Transition

At how small radius can a rolling disk rotate around its precession center? What is the disk's angular velocity and its precessional angular velocity? Let us look at the physics that determines these values.

**At how small radius can a rolling disk rotate around its precession center? What is the disk's angular velocity and its precession rate? Let us look at the physics that determines these values.**

## Given Variables

We assume to know the values of the following variables:

- $R$: The distance from the center of precession to the contact point of the disk
- $r$: The radius of the disk
- $g$: The acceleration due to gravity
- $\theta$: The angle between the disk and the vertical to the ground

## Variables, Units, and Assumptions for Searched Values

We define the following assumptions and variables which we want to calculate:

- $\omega$: The angular velocity (unit rad/s) around the center of mass of the disk
- $\Omega$: The precession rate (unit rad/s), also known as the precession angular velocity, of the disk
- The minimum value of $R$ for the phenomenon to occur
- Lastly, we assume that the disk is homogeneous.

## Relation Between Disk's Angular Velocity and Precession Rate

The speed around the center of mass of the disk, $\omega r$, is the same as that around the precession point, $\Omega R$. This leads to the following:

$$
    \omega=\frac{\Omega R}{a}.
$$

## Equation For Precession Rate

Projecting the force of gravity that acts upon the disk on the local coordinate system shown in the figure above the moment around the contact point becomes the following:

$$
    \vec{M}= r \hat{z} \times mg (-\sin{\theta} \hat{x} - \cos{\theta} \hat{z}) = -mgr \sin{\theta} \hat{y}.
$$

The center of mass is closer to the center of precession than the contact point so the the velocity of the center of mass is

$$
    \vec{v} = \Omega (R - r \sin{\theta}) \hat{y}.
$$

The rotational velocity around the center of mass in the vector in the local coordinate system is

$$
    \vec{\omega} = \begin{bmatrix}
    \frac{\Omega R}{r} \\
    0 \\
    0
    \end{bmatrix}
    .
$$

The precession rate translated to the local coordinate system:

$$
    \vec{\Omega} = \begin{bmatrix}
    \Omega \sin{\theta} \\
    0 \\
    \Omega \cos{\theta}
    \end{bmatrix}
    .
$$

We also need the [moment of inertia around the center of mass of the disk](https://www.engineeringguidance.com/the-moment-of-inertia-of-a-rod-and-a-disk), which is

$$
    I = \begin{bmatrix}
    2mr^2k & 0 & 0 \\
    0 & mr^2k & 0 \\
    0 & 0 & mr^2k
    \end{bmatrix},
$$

where $k=\frac{1}{4}$. We can now use expressions in vector form so the angular momentum around the contact point becomes

$$
    \vec{H} = I (\vec{\omega} + \vec{\Omega}) + r \hat{z} \times \vec{v} = mr \Omega ((r \sin{\theta} - R)(2k +2) \hat{x} + r(k+1) \cos{\theta} \hat{z}).
$$

Then we take the time derivative of the angular momentum at the contact point:

$$
        \frac{d \vec{H}}{dt} = \frac{\partial \vec{H}}{\partial t} + \vec{\Omega} \times \vec{H}
$$

Since the components are not time-dependent, the time-related derivative $\frac{\partial \vec{H}}{\partial t}$ is zero. Therefore the derivative of the angular momentum around the contact point becomes

$$
        \vec{\Omega} \times \vec{H} = -\Omega^2 r m \cos{\theta} (rk \sin{\theta} + (R - r \sin{\theta})(2k + 1)) \hat{y}.
$$

The derivative of the angular momentum at the contact point is equal to the moment at the contact point so

$$
        -mgr \sin{\theta} \hat{y} =  -\Omega^2 r m \cos{\theta} (rk \sin{\theta} + (R - r \sin{\theta})(2k + 1)) \hat{y}.
$$

With that equation, $\Omega$ can be solved and we get the equation for the precession rate:

$$
    \Omega = \sqrt{\frac{g \tan{\theta}}{rk \sin{\theta} + (R - r \sin{\theta})(2k + 1)}}.
$$

## Minimum value of $R$

When the equation above becomes undefined, that is, when the denominator becomes zero, the minimum value of $R$ for precession to occur. So we have have:

$$
    0 = rk \sin{\theta} + (R - r \sin{\theta})(2k + 1).
$$

So if we solve for $R$ we get:

$$
    R = \frac{5}{6} r \sin{\theta}.
$$

## Analysis of Solution

During these calculations, an assumption regarding friction has been made. Some friction must exist to keep the disk in a circular path. Furthermore, these calculations also assume that the disk has no height, which means that $I_{22}$ and $I_{33}$ in $I$ do not have a height term. On top of that, the edge radius in these calculations cannot exist for the length vectors to be correct and at the same time, a larger edge radius would provide a larger contact area, which would lead to greater maximum friction.

If $\theta = 0$, no moment would arise since the gravitational and normal forces would align and precession could not occur. If the velocity were also very low, only a slight instability would cause a moment to arise and since the time derivative of the angular momentum must equal the moment, this means that a large change occurs in the angular momentum when it is small. This would likely lead to the disk falling.

Precession Rate of a Rolling Disk

A charge density is considered spherically symmetric when its value depends only on the distance from a specific point space and not the direction. Consequently, angles do not matter, so it looks the same if the system is rotated. In this post, we will walk through a method of determining 40Ca's nucleus charge density.

**A charge density is considered spherically symmetric when its value depends only on the distance from a specific point space and not the direction. Consequently, angles do not matter, so it looks the same if the system is rotated. In this post, we will walk through a method of determining $^{40}$Ca's nucleus charge density.**

## Charge Density of $^{40}$Ca

As an example, let us determine the spherically symmetric charge density

$$
    \rho_\text{ch}(r;\mathbf{X}) = \frac{\rho^0_{\text{ch}}}{1+e^{(r-b)/a}}
$$

for $^{40}$Ca. Therefore, the constants $\mathbf{X} = [\rho^0_{\text{ch}},a,b]$ needs to be determined. 

## Doubly Magic Nuclei

$^{40}$Ca has a double magic number, meaning the atomic nucleus is more stable than other isotopes. Magic atomic nuclei have very small electric quadrupole moments, making them more or less spherical, which is their most tightly bound state.

## Measuring $^{40}$Ca's Charge Density

This is done by elastically irradiating the target $^{40}$Ca with electrons having an energy of 250 MeV. The target is assumed to be stationary, and no recoil is assumed to occur since the target is much heavier than the electron. Since the target is assumed to be stationary in the LAB frame and no recoil is assumed for the Ca atom, the COM (center-of-momentum) frame does not move in the LAB frame. Thus, energies, scattering cross-sections, and angles are approximately the same in both frames.

Differential cross-sections $\mathcal{E}_i$ has been measured for the electrons to be scattered at an angle $\theta$. These cross-sections have measurement errors $\delta_i$. At the same time, the cross-section can be theoretically expressed as

$$
    \mathcal{T}(\mathbf{X}) = \left(\frac{\sigma}{\Omega} \right)_\text{Mott} |F|^2 = \frac{Z_p^2 Z_t^2 \alpha^2 (\hbar c)^2}{4 \beta^4 E^2 \sin^4(\theta/2)} [1-\beta^2 \sin^2(\theta/2)] \cdot |F|^2,
$$

where $F$ is the form factor, $\alpha$ is the fine structure constant, $\beta = v/c$ with $v$ being the electron's speed, $E$ is the electron's energy, $Z_p$ and $Z_t$ are the charges of the electron and the target, respectively, and $\theta$ is the scattering angle. In order for the measured cross-sections to agree with the theoretical values, minimization was performed on

$$
    \mathcal{X} ^2(\mathbf{X}) = \sum_{i=1}^{N_\theta} \left(\frac{\mathcal{T}(\mathbf{X})-\mathcal{E}_i}{\delta_i} \right)^2.
$$

This needs to be done with the constraint

$$
    Z_t = \frac{4 \pi}{e} \int_0^\infty r^2 \rho_\text{ch}(r;\mathbf{X}) dr.
$$

The constraint ensures that the total charge of the charge density equals the charge of the atomic nucleus $Z_t e$. Without the constraint, the form factor and charge radius would yield incorrect values. The charge density could potentially also become negative, which is not possible. This was done using **scipy**'s **basinhopping**, where a starting point

$$
    \mathbf{X} = [0.06 \text{ e fm}^{-3}, 1.07 \cdot 40^{1/3} \text{ fm}, 0.54 \text{ fm}]
$$

was used.

## Measurements vs Theoretical Value

The vector $\mathbf{X}_*$ that gives the minimum $\mathcal{X}_* ^2$ sum is $\mathbf{X} = [0.077 \text{ e fm}^{-3}, 3.8 \text{ fm}, 0.47 \text{ fm}]$. In the figure below, the charge density can be observed. In the figure below that, a plot of the theoretical and measured cross-sections can also be seen.

![Approximate charge density](https://i.imgur.com/Zr5Vqo0.png "Approximate charge density")

![Cross-section measured vs theoretical](https://i.imgur.com/handYRg.png "Cross-section measured vs theoretical")

## $^{40}$Ca's rms charge radius

Using the charge function, the rms charge radius can be calculated as follows:

$$
    \sqrt{\langle r^2 \rangle} = \sqrt{\frac{4 \pi}{Z_t e} \int_0^\infty r^4 \rho_\text{ch}(r;\mathbf{X}_*) dr} \approx 3.41 \text{ fm}.
$$

This differs from the accepted value of 3.4776(19) fm. The difference may be due to the fact that the charge density function is an approximation. The "true" charge distribution has a similar form, but it does not precisely follow the model used. Furthermore, there is uncertainty in the value of $\rho^0_{\text{ch}}$ at small $r$.

Spherical Symmetric Charge Distribution of Nucleis

The ARMA process, short for autoregressive moving average, is one the most basic time series models but it is an important building block in forecasting. Let us look at its definition and some of its properties such as causality and invertibility but also how to choose which order your ARMA model should have.

**The ARMA process, short for autoregressive moving average, is one the most basic time series models but it is an important building block in forecasting. Let us look at its definition and some of its properties such as causality and invertibility but also how to choose which order your ARMA model should have.**

## ARMA Process Definition

The ARMA($p,q$) process $X_t$ has the following defintion:

$$
    X_t - \sum_{i=1}^q \phi_j X_{t-i} = Z_t + \sum_{i=1}^q \theta_j Z_{t-i},
$$

where $Z_t \sim \text{WN}(0,\sigma^2)$ (white noise with mean 0 and variance $\sigma^2$) and $1 -\sum_{i=1}^q \phi_i z^i$ and $1 + \sum_{i=1}^q \theta_i z^i$ have no common zeros. You might also define the ARMA model in the following short form:

$$
    \phi(B) X_t = \theta(B) Z_t,
$$

where $B$ is the lag operator defined as $B X_t = X_{t-1}$. The two functions are defined as:

$$
    \phi(z) = \sum_{i=1}^q \phi_i z^i
$$

and

$$
    \theta (z) = \sum_{i=1}^q \theta_i z^i
$$

## When is an ARMA Process Stationary?

An ARMA process $X_t$ has a unique (weakly) stationary if and only if $\phi(z) \neq 0$ for all $z \in \mathbb{C}$ when $|z| = 1$

## When is an ARMA Process Causal?

An ARMA process $X_t$ is causal if there exist a real value sequence $(\varepsilon_i, i \in \mathbb{N}_0)$ such that for all $t \in \mathbb{Z}$:

$$
    \sum_{i=0}^\infty |\varepsilon_i| < \infty
$$

and

$$
    X_t = \sum_{i=0}^\infty \varepsilon_i Z_{t-i}.
$$

In other words, $X_t$ can be described as an MA$(\infty)$ process. Equivalently an ARMA process $X_t$ is causal if and only if $\phi(z) \neq 0$ for all $z \in \mathbb{C}$ when $|z| \leq 1$.

## When is an ARMA Process Invertible?

An ARMA process $X_t$ is invertible if there exist a real value sequence $(\pi_i, i \in \mathbb{N}_0)$ such that for all $t \in \mathbb{Z}$:

$$
    \sum_{i=0}^\infty |\pi_i| < \infty
$$

and

$$
    Z_t = \sum_{i=0}^\infty \pi_i X_{t-i}.
$$

In other words, $X_t$ can be described as an AR$(\infty)$ process. Equivalently an ARMA process $X_t$ is invertible if and only if $\theta(z) \neq 0$ for all $z \in \mathbb{C}$ when $|z| \leq 1$.

## Choosing The Order of The ARMA Process

It’s always possible to fit an ARMA$(p,q)$ model with excessively large $p$ and $q$ values, but this isn’t advantageous for forecasting. While this approach often yields a small estimated white noise variance, the mean squared error of forecasts is also affected by errors in parameter estimation. To address this, we introduce a "penalty factor" to discourage the selection of overly complex models.

## AICC Criterion

We start by introducing the AICC criterion, where AIC stands for Akaike’s Information Criterion and the last C denotes biased-corrected. The AICC estimates the Kullback-Leibler divergence, which measures how different one probability distribution is from another, comparing the estimated distribution to the actual data distribution. It assumes that $Z \sim \text{IID} \mathcal{N} (0, \sigma^2)$, but it shows robustness to moderate deviations from normality, such as when $Z_t$ follows a $t$-distribution.

The AICC criterion says that you have to choose $p$, $q$, $\phi_p$ and $\theta_q$ to minimize the following:

$$
    -2 \ln{L(\phi_p, \theta_q, S(\phi_p,\theta_q)/n) + 2n \frac{p+q+1}{n-p-q-2}},
$$

where $\phi_p = (\phi_1, ... , \phi_p)$ and $\theta_p = (\theta_1, ... , \theta_p)$. To apply this criterion in practice, we fit a wide range of models with varying orders (p, q) to the data and select the one that minimizes the negative log-likelihood, adjusted by the penalty factor $2n \frac{p+q+1}{n-p-q-2}$

## BIC Criterion

One issue with the AICC criterion is that the estimators for p and q aren’t consistent; they don’t converge almost surely to the true values. In contrast, consistent estimators can be derived using the Bayesian Information Criterion (BIC), which also penalizes the selection of large p and q values, helping to prevent overfitting.

The BIC criterion says that you have to choose $p$ and $q$ to minimize the following:

$$
    (n-p-q) \ln{\frac{n\hat{\sigma}^2}{n-p-q}}+n (1+\ln{\sqrt{2 \pi}}) + (p+q) \ln{\frac{\sum_{t=1}^n X_t^2 - n \hat{\sigma}^2}{p-q}},
$$

where $\hat{\sigma}^2$ is the maximum likelihood estimate of the white noise variance. A downside of the BIC is its efficiency in finding minimizers. While selecting the model that minimizes the AICC is asymptotically efficient for causal and invertible ARMA processes, this isn’t the case for the BIC. In this context, efficiency means that minimizing the AICC will, in the long run, lead to a model with the lowest one-step-ahead prediction errors.

When using constrained maximum likelihood estimation—where certain coefficients are assumed to be zero during the estimation process—the term $p+q+1$ is replaced by $m$, which represents the number of non-zero coefficients, in both the AICC and BIC.

The ARMA Process

S-parameters, short for scattering parameters, are a useful and easy way to measure a network's response to input at its terminals (ports). With the S-parameters of a network determined, the behavior of the network in a new external environment can be determined. What is measured is the reflected and transmitted waves going through the network's ports.

**S-parameters, short for scattering parameters, are a useful and easy way to measure a network's response to input at its terminals (ports). With the S-parameters of a network determined, the behavior of the network in a new external environment can be determined. What is measured is the reflected and transmitted waves going through the network's ports.**

## S-parameter Formula

To calculate the S-parameters of a network one can use the normalized incident complex voltage wave $a_i$ at the $i$:th port defined as

$$
    a_i = \frac{V_i + Z_i I_i}{2 \sqrt{|\text{Re} Z_i|}}
$$

and the normalized reflected complex voltage wave $b_i$ defined as

$$
    b_i = \frac{V_i - \overline{Z_i} I_i}{2 \sqrt{|\text{Re} Z_i|}},
$$

where $V_i$ is the terminal voltage, $I_i$ the terminal current and $Z_i$ is an arbitrary reference impedance (and $\overline{Z_i}$ is it's complex conjugate). Then the S-parameter $S_{ij}$ is defined as

$$
    S_{ij} = \frac{b_i}{a_j}.
$$

## S-parameters For a 2-port Network

For a 2-port network, seen in the figure above, there are four different S-parameters. These can be used to describe the relation between incident waves $a_1$ and $a_2$ and the reflected waves $b_1$ and $b_2$ in the following

$$
    \left\{
    \begin{array}{cc}
         b_1 = & S_{11} a_1 + S_{12} a_2 \\
         b_2 = & S_{21} a_1 + S_{22} a_2
    \end{array}
    \right.
    .
$$

To get $S_{11}$ the second port has to be terminated with a matched load $Z_L = Z_0$ and $a_2 = 0$ so

$$
    S_{11} = \left. \frac{b_1}{a_1} \right|_{a_2 = 0}.
$$

To get $S_{22}$ the first port has to be terminated with a matched load $Z_S = Z_0$ and $V_S = 0$ so

$$
    S_{22} = \left. \frac{b_2}{a_2} \right|_{a_1 = 0}.
$$

To get $S_{21}$ the second port has to be terminated with a matched load $Z_L = Z_0$ and $a_2 = 0$ so

$$
    S_{21} = \left. \frac{b_2}{a_1} \right|_{a_2 = 0}.
$$

To get $S_{12}$ the first port has to be terminated with a matched load $Z_S = Z_0$ and $a_1 = 0$ so

$$
    S_{12} = \left. \frac{b_1}{a_2} \right|_{a_1 = 0}.
$$

Using these formulas and by choosing a real reference impedance $Z_0$ (so $Z_0=\overline{Z_0}$ we get

$$
    S_{11} = \frac{b_1}{a_1} = \frac{\frac{V_1}{I_1} - Z_0}{\frac{V_1}{I_1} + Z_0} = \frac{Z_1 - Z_0}{Z_1 + Z_0}.
$$

Therefore we also get

$$
    Z_1 = Z_0 \frac{1 + S_{11}}{1 - S_{11}},
$$

where at port 1, $Z_1 = \frac{V_1}{I_1}$ is the input impedance.

## S-parameters For Antennas

![S-parameters For Antennas](https://i.imgur.com/sRWuN8i.png "S-parameters For Antennas")

S-parameters for an $N$:th element antenna array is calculated by

$$
    S_{mn} = \left. \frac{V_m^-}{V_n^+} \right|_{V_k^+ = 0 \text{ for } k \neq n}
$$

where $V_m^+$ are the incident and $V_m^-$ the reflected voltage wave amplitudes for the $m$:th array element seen in the figure above.

S-parameters in Circuits and Antennas

In this post, vehicle dynamics parameters are described which can affect the aerodynamic performance of a car.

**In this post, vehicle dynamics parameters are described which can affect the aerodynamic performance of a car.**

## Steering Angle

The steering angle is the angle between the direction of a wheel and the car's longitudinal axis. The steering angle for each front wheel can be different if the car has Ackermann steering or anti-Ackermann steering. An illustration can be found in the figure above.

## Body Slip Angle

![Body slip angle](https://i.imgur.com/rZFkoDP.png "Body slip angle")

The body slip angle is a yaw angle meaning it is the angle that is created when the car rotates around its vertical axis. The rotational axis goes through the car's center of gravity. The body slip angle is the angle between the velocity vector of the car and the car's longitudinal direction. A visual representation of the body slip angle can be seen in the figure above.

## Pitch Angle

![Pitch angle](https://i.imgur.com/gckbgDJ.png "Pitch angle")

The pitch angle is the angle created if the vehicle's sprung mass (everything but the wheels) rotates around its lateral axis, as shown in the figure above. The vehicle rotates around what is called its pitch center.

## Roll Angle

![Roll angle](https://i.imgur.com/RsVYKF6.png "Roll angle")

The roll angle is the angle that is created if the vehicle's sprung mass rotates around its longitudinal axis, see the figure above. The vehicle rotates around what is called its roll center.

## Ride Height

![Ride height](https://i.imgur.com/iGjAZla.png "Ride height")

The ride height is the smallest distance between the ground and the car's floor when it is stationary, see the figure above.
    
## Aerodynamic Influence on Vehicle Dynamics

Vehicle dynamics is a large and complex subject. Mass distribution, tires, suspension setup, and aerodynamics are examples of areas that greatly affect the dynamics of a car.

One of the main aims of aerodynamics on a race car is to increase its negative lift force (downforce). Greater downforce makes it possible for the tires to produce larger forces (grip) while accelerating/braking (longitudinal forces) and cornering (lateral forces). The basic formula describing this relationship is $F_{\text{long/lat}} = F_z \mu$, where $F_{\text{long/lat}}$ is the longitudinal or latitudinal force, $F_z$ the normal force which is increased with increased downforce and $\mu$ the friction coefficient. Consequently, the speed of the car can be increased in most driving scenarios or the same grip can be achieved with less tire slip, i.e. smaller slip angle, when compared to a car with less downforce. This helps reduce heat production and therefore preserve tire life. Furthermore, smaller slip angles are easier to control for the driver. The car's stability in this case is improved compared to larger slip angles and aerodynamic parts tend to work more efficiently, improving the balance.

In most cases, a weight and aerodynamic distribution close to 50/50 on the car's rear and front axle is desired. A car with the center of gravity (COG) closer to the front axle usually generates more understeer. Conversely, the aerodynamic center of pressure (COP or aero balance), the point at which the resultant downforce is located, closer to the front usually generates more oversteer. A car with a front-biased COG can therefore be counterbalanced by placing the COP forward. However, it is favorable in some cases to design the aerodynamic package with a COP slightly rear-biased. This can generate understeer, but understeer is considered more stable and easier to control when compared to oversteer.

## The Front Wing

Vehicle dynamics parameters described above are closely related to the aerodynamic package on a given car. While the generated downforce affects body slip, pitch, and roll angles, these parameters in turn affect the aerodynamic efficiency. Pitch angles for example can have a large impact on vehicle balance. When braking and pitching forwards, the front wing's/front splitters's downforce is generally increased. 

However, other factors such as already low static ride height or an extreme pitch angle might cause a decrease in downforce for the front wing. This also results in less airflow to the car's floor and diffuser, decreasing the middle and rear downforce. If the front wing downforce increases, the aerodynamic balance shifts forward, decreasing the stability of the car. While this is considered unfavorable, increased grip on the front axle also results in improved braking capabilities. On the contrary, if the downforce on the front wing decreases the balance shifts rearwards and the opposite occurs. It is more difficult to describe the influence of steering and roll angles. These affect the aerodynamics in an asymmetric manner and can cause unpredictable behavior in the car. The same is true for body slip angles, or yaw angles, where the airflow hits the car with an angle from the longitudinal axis.

## The Diffuser

A diffuser expands the airflow underneath the car's floor. This expansion forces the flow upstream to increase in speed in order to fill up the larger volume in the diffuser, decreasing the pressure upstream. In motorsport, it has been observed that diffusers are particularly effective aero devices. The reason is that diffusers increase the vehicle's downforce under the whole floor. Downforce in this region is close to the center of the car and distributed over a large area. This results in a good balance contribution and makes the vehicle less sensitive to pitch or other phenomena that affect the balance. 

In order for a diffuser to work efficiently, the front wing needs to supply the floor and diffuser area with sufficient airflow. Optimally the airflow's velocity is large and laminar. Another option for increasing the airflow is by feeding the diffuser with air from the sides of the vehicle. This can however interfere with the low-pressure diffuser region and increase the pressure, resulting in a decrease of downforce. A compromise between keeping the floor sealed and increasing the airflow is usually optimal. This airflow can also be used to induce vortices.

Relation Between Vehicle Dynamics and Aerodynamics

Signal distortion is unwanted changes in the output signal relative to the input signal. Some of these distortions can be harmonic, meaning that disturbances occur at multiples of the fundamental frequency.

**Signal distortion is unwanted changes in the output signal relative to the input signal. Some of these distortions can be harmonic, meaning that disturbances occur at multiples of the fundamental frequency.**

## Total Harmonic Distortion

Distortion can be measured in various ways, but a common measurement is total harmonic distortion (THD). THD is defined as

$$
    \text{THD} = \frac{\sqrt{V_2^2+V_3^2+ ... +V_n^2}}{V_1},
$$

where $V_n$ is the n-th harmonic amplitude in the frequency domain. $V_1$ is called the fundamental frequency. Generally, only the first five or six harmonic amplitudes are significant for THD, and often the first three are sufficient in practice without greatly impacting the margin of error.

## Slew Rate of Operational Amplifier

Harmonic distortion manifests in various forms. One example is the slew rate of an operational amplifier. The slew rate denotes the maximum rate of change of the output signal, dictated by the specified slew rate value of the operational amplifier. When this rate of change is exceeded, the waveform is replaced by a linear ramp that adheres to the maximum slew rate until it eventually intersects with the original waveform again, as illustrated in the figure above. This phenomenon can lead to significant alterations in signal fidelity, particularly in high-frequency applications, where the ability to accurately track rapid changes in the input signal is critical for maintaining overall performance.

## Clipping in Operational Amplifier

Another form of harmonic distortion is clipping. Clipping can be categorized into soft clipping and hard clipping. 

### Hard Clipping

![Hard Clipping](https://i.imgur.com/GlgkIKr.png "Hard Clipping")

Hard clipping means that the output signal can only reach a certain amplitude, resulting in the waveform appearing "chopped off." For a negatively feedbacked operational amplifier, this amplitude depends, among other factors, on the power supply voltage. The output voltage cannot exceed the operational amplifier's supply voltage. Another way to achieve hard clipping is by using diodes between the output voltage and ground, as the output voltage in that case cannot exceed the threshold voltage

### Soft Clipping

![Soft Clipping](https://i.imgur.com/Bg2ygds.png "Soft Clipping")

Soft clipping, on the other hand, means that the peaks of a waveform experience a reduced increase, resulting in a more "rounded" appearance. Soft clipping in an op-amp occurs when the output signal approaches the supply voltage limits but does not abruptly cut off; instead, it gradually compresses the signal as it nears those limits. This typically happens when the input signal exceeds the linear operating range of the op-amp, causing the output to start saturating.

Harmonic Distortion in Operational Amplifiers

The Angstrom Method is a method to determine the thermal conductivity or the specific thermal conductivity of different materials. The method builds upon sending temperature pulses from a reservoir and measuring the temperature at different points in the material.

**The Angstrom Method is a method to determine the thermal conductivity or the thermal diffusivity of different materials. The method builds upon sending temperature pulses from a reservoir and measuring the temperature at different points in the material.**

## The Heat Equation

The situation is determined by the heat equation,

$$
    \frac{\partial T}{\partial t}(\mathbf{r},t) = \alpha \Delta T(\mathbf{r},t) - \epsilon T(\mathbf{r},t),
$$

where $T(\mathbf{r},t)$ is the temperature as a function of position $\mathbf{r}$ and time $t$, $\alpha$ is the thermal diffusivity of the material and $-\epsilon T$ represents the heat losses to the ambient air. In the case of the material being an oblong rod, the situation can be approximated as one-dimensional, which gives the following formulation of the heat equation

$$
    \frac{\partial T}{\partial t}(x,t) = \alpha \frac{\partial^2 T}{\partial x^2}(x,t) - \epsilon T(x,t).
$$

Furthermore, $\alpha = \frac{\kappa}{c \rho}$, where $\kappa$ is the thermal conductivity, $s$ is the speific heat capacity and $\rho$ is the density, as well as $\epsilon = \frac{R O}{S c \rho}$ where $R$ is an emission constant, $O$ is the circumference of the rod and $S$ the cross-sectional area of the rod, meaning

$$
    \frac{\partial T}{\partial t}(x,t) = \frac{\kappa}{c \rho} \frac{\partial^2 T}{\partial x^2}(x,t) - \frac{R O}{S c \rho} T(x,t)
$$

According to the Angstrom method, a heater shall now be connected to one of the rod's ends, say at $x=0$, at which it works as a reservoir securing the boundary condition $T(0,t) = T_0 \cos \omega t$ for any chosen temperature $T_0$ och ''angular velocity'' $\omega$. Due to diffusion and emission, the temperature $T$ will not only vary sinusoidal with $t$ and $x$ but also exponentally decrease with $x$. Therefore, a solution with the form

$$
    T(x,t) = A e^{-ax} \cos(\omega t - b x),
$$

## Using Two Measuring Points to Determine $\kappa$

where $A$, $a$ (probably $>0$) and $b$ are constants, will be assumed. If this is put into the heat equation for the oblong rod then

$$
    -\omega A e^{-a x} \sin(\omega t - b x) = \frac{\kappa}{c \rho} A e^{-a x} \big((a^2-b^2) \cos(\omega t - b x) - 2 a b \sin(\omega t - b x)\big) - \frac{R O}{S c \rho} A e^{-a x} \cos(\omega t - b x),
$$

is required, which if the exponential factor $A e^{-a x}$ is canceled, gives two equations in $\sin$ and in $\cos$

$$
    \sin \implies -\omega = 2 a b \frac{\kappa}{c \rho},
$$

$$
    \cos \implies \frac{\kappa}{c \rho} (a^2-b^2) - \frac{R O}{S c \rho} = 0.
$$

Furthermore, the Angstrom method is based on measuring the temperature, $T_1$ and $T_2$, at two points, $x_1$ and $x_2$ (with $x_2>x_1$ for simplicity), and by comparing the amplitudes, say $A_1$ and $A_2$ respectively, at those points. With the nature of the temperature function, it follows that

$$
    \frac{T_2}{T_1} = \frac{A_2}{A_1} = e^{-a(x_2-x_1)},
$$

or 

$$
    a = \frac{\log(\frac{A_1}{A_2})}{x_2-x_1},
$$

and by comparing the phases, say $\varphi_1$ and $\varphi_2$ respectivaly, it follows that

$$
    \varphi_2-\varphi_1 = b(x_2-x_1),
$$

or

$$
    b = \frac{\varphi_2-\varphi_1}{x_2-x_1}.
$$

and with these equations, $a$ and $b$ can be eliminated, which gives the following expression for the thermal conductivity $\kappa$

$$
    \kappa = \frac{c \rho \omega (x_2-x_1)^2}{2 (\varphi_1-\varphi_2) \log(\frac{A_1}{A_2})}.
$$

## Determining Thermal Diffusivity

Using the above equation for $\kappa$ we can insert that into the equation for the thermal diffusivity

$$
    \alpha = \frac{\kappa}{c \rho} = \frac{\omega (x_2-x_1)^2}{2 (\varphi_1-\varphi_2) \log(\frac{A_1}{A_2})}.
$$

So that is how you determine the thermal diffusivity and thermal conductivity of a material using the Angstrom Method.

Determining Thermal Conductivity Using The Angstrom Method

To apply models like the ARMA model to time series sample data, the seasonality and trend have to be removed. Let's look at how to do that.

**To apply models like the ARMA model to time series sample data, the seasonality and trend have to be removed. Let's look at how to do that.**

Let us say we have time series $X$ which has both a trend and a seasonality then we can do a decomposition of the time series and write it as:

$$
    X_t = m_t + s_t + Y_t,
$$

where $m_t$ is the trend component, $s_t$ is the seasonal component with period $d$ and $Y_t$ is the rest which is [weakly stationary](https://www.engineeringguidance.com/time-series-white-noise) and has mean of 0. Since we have singled out a seasonal component we know that $s_t = s_{t+d}$ because $d$ is the period of the seasonality.

Before we discuss how to eliminate the trend and seasonal component, we have to examine a couple of time series operators. 

## The Lag Operator

The lag operator $B$, also called the backshift operator, is defined as:

$$
    B X_t = X_{t-1}.
$$

So when taking the lag operator to the power of $k$ we get:

$$
    B^k X_t = X_{t-k}.
$$

## The Difference Operator

The difference operator $\nabla$ is defined as:

$$
    \nabla X_t = (1-B) X_t =  X_t - X_{t-1}.
$$

So when taking the difference operator to the power of $k$ and using the binomial theorem we get:

$$
    \nabla^k X_t = (1-B)^k X_t = \sum_{i=0}^k {k \choose i} (-1)^i B^i X_t = \sum_{i=0}^k {k \choose i} (-1)^i X_{t-i}.
$$

## The Lag-d Differencing Operator

The lag-d differencing operator $\nabla_d$ is defined as:

$$
    \nabla_d X_t = (1-B^d) X_t =  X_t - X_{t-d}.
$$

## Removing Seasonality in a Time Series

To eliminate seasonality in a time series we will have to use the lag-d differencing operator. Since the seasonality has a period of $d$ then:

$$
    \nabla_d X = m_t - m_{t-d} + s_t - s_{t-d} + Y_t - Y_{t-d} =  \nabla_d m_t + \nabla_d Y_t,
$$

because $s_t = s_{t-d}$.

## Removing Trend in a Time Series

Eliminating trend in a time series is a bit more tricky than removing seasonality. We have to assume that

$$
    m_t = \sum_{i=0}^q a_i t^i,
$$

where $a_i$ is some real constant and $q < n$ where $n$ is the sample size of the data. Lets apply the difference operator to the power of $q$ on $m_t$:

$$
    \nabla^q m_t = \sum_{i=0}^q {q \choose i} (-1)^i \sum_{j=0}^q a_j (t-i)^j.
$$

Since the sums are finite we can change the order of them, so we get the following

$$
    \sum_{j=0}^q a_j \sum_{i=0}^q {q \choose i} (-1)^i (t-i)^j.
$$

Now, we focus on the inner sum:

$$
    \sum_{i=0}^q {q \choose i} (-1)^i (t-i)^j.
$$

By applying the binomial theorem and properties of polynomial expansions on $(t-i)^j$ we get:

$$
    (t-i)^j = \sum_{k=0}^j {j \choose k} t^{j-k} (-i)^k.
$$

By substituting this into our sum and changing order of the sums we get:

$$
    \sum_{i=0}^q {q \choose i} (-1)^i (t-i)^j = \sum_{k=0}^j {j \choose k} t^{j-k} (-1)^k \sum_{i=0}^q {q \choose i} (-1)^i i^k.
$$

Using the fact that

$$
    \sum_{i=0}^q {q \choose i} (-1)^i i^k = 0 \text{ for } k < q
$$

we find that the entire expression evaluates to 0 when $j < q$. When $j = q$ this inner sum evaluates to $q!$. Therefore, the final expression simplifies to:

$$
    a_q q!
$$

and this means that

$$
    \nabla^q m_t = a_q q!.
$$

So when we apply $\nabla^q$ to $X_t$ we get

$$
    \nabla^q X_t = a_q q! + \nabla^q Y_t.
$$

Given that $Y$ is regarded as a stationary process with a mean of zero, it can be demonstrated that this property extends to $\nabla^q Y$. Consequently, if $s = 0$, $\nabla^q X$ behaves as a stationary process with a mean of $a_q q!$.

Time Series Seasonality and Trend

The patch antenna is the component used to connect the Transmit Receive Module's amplifiers to free space, allowing for electromagnetic waves to be transferred with minimal losses. For most antennas to work efficiently, however, the radiated electromagnetic energy not only needs to be transmitted but also directed in order to be able to transmit efficiently and over longer distances.

**The patch antenna is the component used to connect the [Transmit Receive Module's](https://www.engineeringguidance.com//what-is-a-transmit-receive-module) amplifiers to free space, allowing for electromagnetic waves to be transferred with minimal losses. For most antennas to work efficiently, however, the radiated electromagnetic energy not only needs to be transmitted but also directed in order to be able to transmit efficiently and over longer distances.**

## What is an Antenna Array?

A single antenna element usually has a wide radiation pattern but for many applications such as for radars, it is necessary to have antennas with a high directivity that can create beams that focus the radiated power. To do this the most common way is to form an array of antenna elements, known as an antenna array. The total radiated field of such an array is determined by how the individual elements' radiation interfere in different directions which results in constructive interference in some directions while others experience destructive interference. 

## Two Element Antenna Array

![Two element linear array](https://i.imgur.com/CruAuUo.png "Two element linear array")

The smallest antenna array possible is a simple two element linear array as the one shown in the figure above consisting of two infinitesimal horizontal dipoles positioned at $\pm d/2$ along the z-axis. Assuming no coupling between the elements the radiated field is given by adding the field of the two individual antenna elements giving

$$
    \mathbf{E}_t = \mathbf{E_1} + \mathbf{E_2} = \hat{\theta}j\eta\frac{kI_0l}{4\pi}\left( \frac{e^{-j[kr_1-(\beta/2)]}}{r_1}\cos{\theta_1} + \frac{e^{-j[kr_2+(\beta/2)]}}{r_2}\cos{\theta_2}\right).
$$

where $\eta$ is the free space impedance, $k$ is the wavenumber, $I_0$ is the current amplitude, $l$ is the length of the dipoles, and $\beta$ is the phase of the signal. With the assumptions that $\theta\simeq\theta_1\simeq\theta_2$, the phase variations $r_1\simeq r-\frac{d}{2}\cos{\theta}$, $r_2\simeq r+\frac{d}{2}\cos{\theta}$ and the amplitude variation  $r_1\simeq r_2 \simeq r$ the equation reduces to 

$$
    \mathbf{E_t} = \hat{\theta}j\eta\frac{kI_0le^{-jkr}}{4\pi r}\cos{\theta}\cdot\left\{ 2\cos\left[\frac{1}{2}(kd\cos{\theta}+\beta)\right] \right\}.
$$

The last factor within curly brackets is called the array factor (AF). It can be seen by comparing the first and second equations that the radiated field from the array is the product of the radiated field from a single element multiplied by the AF. By varying the separation $d$ and/or the phase $\beta$ the total field of the array can be controlled. Most antennas however have a fixed separation by design and what is left to control is the phase. When controlling the phase for the antenna array it is then called a phased antenna array. 


![Beam angle antenna](https://i.imgur.com/zw4w0mC.png "Beam angle antenna")

The same principle as for the two element array still applies for arrays when extending it to higher dimensions and adding more elements. Using more elements can get even better directivity allowing for beams that can reach longer. An illustration of the beam angles $(\theta,\varphi)$ in three dimensions can be seen in the figure above where it can be seen that at $\theta = 0$ and $\varphi = 0$ the antenna array is sending radiation straight ahead. To get more accurate results, however, the interaction between the elements, known as mutual coupling, has to be accounted for.

## Mutual Coupling

Mutual coupling in an antenna array is a term that is used to describe the interaction between antenna elements. It occurs mainly as a consequence of nearby elements receiving some of the radiated energy but can also occur because of reflections, reradiation, and scattering of electromagnetic radiation. This phenomenon exists for all types of arrays but becomes even more important to account for when working with large arrays. What happens is that when an element is excited by a current, the surrounding antenna elements will receive some of the energy leading to changes in their own respective excitation. This sums up to a total change in the excitation for all elements which will change the characteristics of the antenna seen from the driving source. One key factor is that the antenna's terminal impedance, as seen by the source, changes which might cause an impedance mismatch resulting in a decrease in radiated energy. While designing an antenna array the coupling is accounted for and the antenna array is typically impedance matched for a beam angle straight ahead. To assess the beamforming capability, however, all angles of interest have to be evaluated since the coupling tends to be worse for large beam angles.

## The Patch Antenna

A patch antenna is a type of microstrip antenna that is lightweight, easy to manufacture, and commonly used for a variety of applications. It typically consists of three layers; a ground plane, a substrate, and the patch elements. A patch antenna designed on a printed circuit board (PCB) requires a large grounded plane below the antenna elements with the substrate in between. The signal can be fed to the elements in different ways such as through a microstrip feed line or with a coaxial cable with its feed pin going through the substrate and soldered to the patch, the latter being shown in the figure above.

When designing a patch antenna array, lots of different aspects have to be taken into account. For example, the shape of the patch, the type of material, and the dimensions all have significant effects on the antenna array performance. 

## Design of Patch Antenna Array

While the layout of the patch can take different shapes to better suit a specific application, a rectangular patch antenna with probe feeding is shown in the figure above because of its simplicity. This shape has also proven to have better overall characteristics compared to other simple shapes such as triangular and circular patches.

For a large antenna array, most of the antenna elements will behave similarly with the exception of the elements close to the array edge. Therefore a practical way to design an antenna array is to use a unit cell representing one element within an infinite array. Using a unit cell reduces the simulation time and computational complexity significantly. It is a decomposition-based finite element technique for modeling repetitive antenna arrays and can be used because of the similarities of the electromagnetic fields for elements inside a large array. 

![Hex cell patch antenna](https://i.imgur.com/hgFSM7X.png "Hex cell patch antenna")

The figure above shows a hexagonal unit cell of one of the patch antennas in an infinite array. A hexagonal unit cell can fit more elements within the same amount of area compared to a square cell. The element spacing in all three directions is chosen to be half a wavelength ($\lambda/2$). Using $\lambda/2$ spacing of the elements is chosen since it gives a good trade off between directivity and stearability with low risks of grating lobes. Grating lobes are beams that can be created for larger scan angles if the spacing is higher which would cause a decrease in transmitted power in the desired direction.

Depending on the frequency at which the antenna is designed the remaining parameters such as cell width, patch width, feed pin offset, feed pin radius, and the feed pin's teflon thickness for the antenna element can optimized for this single frequency using a finite element analysis software and common measurements from manufacturers.

Patch Antenna Arrays

In fluid dynamics, there are many concepts, formulas, and coefficients. In this post the most fundamental ones are explained

**In fluid dynamics, there are many concepts, formulas, and coefficients. In this post the most fundamental ones are explained**

## Bernoulli's Principle

In incompressible isentropic flow, i.e. a flow where a negligible amount of energy is lost to non-reversible processes such as heat in turbulent flows, one can derive Bernoulli's principle from the conservation of energy

$$
    \frac{u^2}{2}+\frac{p}{\rho}+z\cdot g= \text{constant},
$$

where $z$ is the vertical position of a fluid element, this law states that an increase in velocity implies a decrease in pressure and vice versa. 

## Dimensionless parameters
Dimensionless parameters are used to make a parameter comparable with parameters of other types of flow. Below the relevant parameters are presented.

### Reynolds Number
The Reynolds number, Re, is a dimensionless number that tells the relation between the inertia and viscosity in a newtonian fluid and is defined by the following equation

$$
    \text{Re} = \frac{\rho V L}{\mu} = \left\{ \frac{\text{Inertia}}{\text{Viscosity}} \right\}.
$$

$V$ and $L$ are the characteristic velocity and length of the flow. A high Re is associated with a fast, large-scale turbulent flow, whereas a low Re corresponds to a slow, viscous flow. Fast flows of gas imply a relatively high Re.

### Drag Coefficient, $C_D$

The drag coefficient of a body, $C_D$, is a dimensionless number indicative of drag force. $C_D$ is defined by the following equation

$$
    C_D = \frac{F_D}{\frac{1}{2}\rho u_{\infty}^2 A},
$$

where $A$ is the projected area from the front of the body and the denominator $(\frac{1}{2} \rho  u_{\infty}^{2})$ is the dynamic pressure of the free stream.

### Lift Coefficient, $C_L$

The lift coefficient of a body, $C_L$ is defined by the following equation

$$
    C_L = \frac{F_L}{\frac{1}{2}\rho u_{\infty}^2 A},
$$

where $A$ is the projected area from the front of the body. $C_L$ is a dimensionless coefficient that increases with lift force. Depending on the area of study, $C_L$ will have different directions. For example, when designing race cars, $C_L$ will indicate the lift in the negative $ z$ direction, i.e., a high downforce will correspond to a large $C_L$. When designing airplanes it will be in the other direction to indicate lift.

### Pressure Coefficient, $C_p$

The pressure coefficient, $C_p$, is a dimensionless parameter that describes the relative pressure. $C_p$ is defined as 

$$
    C_{p}=\frac{p-p_{\infty}}{\frac{1}{2} \rho  u_{\infty}^{2}},
$$

where $p$ is the static pressure at the point where $C_p$ is calculated, while $p_{\infty}$ is the static pressure in the free stream. 

### Skin Friction Coefficient, $C_f$

The skin friction coefficient, $C_f$, is a dimensionless parameter defined as

$$
    C_{f}=\frac{\tau _{w}}{\frac{1}{2} \rho  u_{\infty}^{2}}.
$$

### Vorticity

The vorticity, $\vec{\zeta}$, is defined as the curl of the velocity field,  $\nabla\,\times\,\vec{u}$, which relates to the rotation of a velocity field. If the vorticity is zero the flow has no rotation and is called irrotational. The vorticity in the $x$-direction is given by 

$$
    \zeta_x= \frac{\partial w}{\partial y} - \frac{\partial v}{\partial z}.
$$

## Navier-Stokes Equations

A fluid can be described with the Navier-Stokes equations if it is incompressible and newtonian as 

$$
    \rho \left( \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla) \vec{u} \right) = \rho \vec{g} - \nabla p + \mu \Delta \vec{u}. 
$$

### RANS

In turbulent flow, some fluctuations cause rapid changes in the velocity and pressure in the Navier-stokes equations. These rapid changes can be easier to account for if the velocity and pressure are split into a mean variable $\bar{u}$ and a fluctuation variable $u'$. The velocity in the $x$-direction is therefore $u = \bar{u} + u'$. If those splits are put into the Navier-stokes equations we get the Reynold's Averaged Navier Stokes (RANS) equations after some derivations. The RANS in $x$-direction is 

$$
     \rho\frac{d \overline{u}}{dt}=-\frac{\partial \overline{p}}{\partial x}+\rho g_x+\frac{\partial }{\partial x}\left(\mu\frac{\partial \overline{u}}{\partial x}-{\rho \overline{u^{\prime 2}}}\right)+\frac{\partial }{\partial y}\left(\mu\frac{\partial \overline{u}}{\partial y}-{\rho \overline{u^{\prime}v^{\prime}}}\right)+\frac{\partial }{\partial z}\left(\mu\frac{\partial \overline{u}}{\partial z}-{\rho \overline{u^{\prime}w^{\prime}}}\right).
$$
 
## Flow Past Boundary

Bodies, such as plates, immersed in a fluid stream create a boundary layer flow. The boundary layer is defined as the region where the flow velocity is less than 99 percent of the velocity of the external flow ($u_\infty$). The boundary layer flow is made up of two parts, the laminar part and the turbulent part:

$$
    \frac{\delta}{x} \approx
    \left\{
    \begin{aligned}
        & \frac{5}{\text{Re}_x^{1/2}} \qquad 10^3 < \text{Re}_x < 10^6 \text{ laminar flow} \\
        & \frac{0.16}{\text{Re}_x^{1/7}} \qquad 10^6 < \text{Re}_x \text{ turbulent flow}
    \end{aligned}
    \right. .
$$

$\delta$ is the thickness of the boundary layer and $\text{Re}_x$ the Reynolds number at $x$. An illustration of he boundary level flow can be seen in figure below.


![Boundary layer flow](https://i.imgur.com/HpCZOT0.png "Boundary layer flow")

### Flow Separation

Backflow, i.e. flow going in the opposite direction of the freestream, at the wall indicates that the flow is separated. The point of separation is when the wall shear is 0 which means the velocity gradient is 0.

### The Logarithmic Overlap Law

The turbulent flow near a wall can be broken up into three regions based on which shear stress dominates:

* Wall layer: Viscous shear dominates
* Outer layer: Turbulent shear dominates
* Overlap layer: A mix of the shear types

The inner layers velocity can be modeled as $u = f(\mu, \tau_w, \rho, y)$, $\tau_w$ is the wall shear stress. By dimensional analysis the following relation can be deducted:

$$
    u^+ = \frac{u}{u^*} = F \left(\frac{y u^*}{\nu} \right),
$$

where $u^* = \sqrt{\tau_w/\rho}$. An akin method can be done with the outer layer, $(u_\infty - u)_{\text{outer}} = g(\delta, \tau_w, \rho, y)$ and the dimensional analysis leads to:

$$
    \frac{u_\infty - u}{u^*} = G\left( \frac{y}{\delta} \right).
$$

The overlap layer makes it so these two layers overlap smoothly and for this to be true the overlap velocity must vary logarithmically with y:

$$
    \frac{u}{u^*} = \frac{1}{\kappa} \ln{\frac{y u^*}{\nu}} + B,
$$

where $\kappa \approx 0.41$ and $B \approx 5.0$. 
We define the dimensionless number $y^+$ as

$$
    y^+ =\frac{y u^*}{\nu},
$$

where $u^* = \sqrt{\tau_w/\rho}$. For $y^+ \in [0,5]$ the inner layer dominates and is proportional to $y^+$. It can be described with the dimensionless variable

$$
    u^+ = \frac{u}{u^*} = y^+.
$$

For $y^+ \in [30, 10^3]$ the overlap layers dominates and can be described as

$$
    u^+ = \frac{1}{\kappa} \ln{y^+} + B,
$$

where $\kappa \approx 0.41$ and $B \approx 5.0$. For $y^+ \in [5,30]$ the inner layer has to curve so it merges with the overlap layer. Only experimental data exists for this interval making it difficult to model.

Fundamental Fluid Dynamics Theory

What is a varactor? The varactor, which stands for variable capacitor, is a diode with a well defined capacitance that can be controlled with the reverse voltage applied to it.

**What is a varactor? The varactor, which stands for variable capacitor, is a diode with a well defined capacitance that can be controlled with the reverse voltage applied to it.**

## How does a varactor diode work?

The varactor's depletion region width $d$, seen in the illustration of the varactor diode in the figure above, increases when the reverse voltage over the diode increases and can be likened to the distance between the plates of a plate capacitor being increased as seen in 

$$
    C = \frac{A \varepsilon}{d},
$$

where $C$ is the capacitance, $A$ is the cross section area of the plate capacitor, $\varepsilon$ the dielectric constant and $d$ is the distance between the two plates. The increase of $d$ due to the increase of reverse voltage causes the capacitance of the varactor to decrease. These characteristics enable the varactor's capacitance to be controlled and set to a specific value. All varactor diodes do however have a maximum allowed reverse voltage known as the breakdown voltage. Exceeding the breakdown voltage could break the varactor which is why the reverse voltage must be kept below this limit.

## Varactor symbol

![Varactor Diode Symbol](https://i.imgur.com/sTzVizV.png "Varactor Diode Symbol")

The varactor's circuit can be seen in the figure above. As can be seen, it is a combination of the symbol of a diode and a capacitor.

## Varactor Modeling in Circuits

The figure below shows the SPICE (Simulation Program with Integrated Circuit Emphasis) model used for modeling the varactor as an equivalent circuit. It is derived from *Skyworks Solutions Inc*'s datasheets. The SPICE model contains parasitic components $C_p$ and $R_s$ and the package inductance $L_s$. The difference compared to *Skyworks Solutions Inc*'s datasheet is the change from the diode model to a simple variable capacitance $C_\text{var}$.

![Linear SPICE model of varactor diode](https://i.imgur.com/Yj6lro4.png "Linear SPICE model of varactor diode")

## Nonlinear Varactor Modeling

While the SPICE model shown in the figure above works for simulating S-parameters, a nonlinear model that accounts for the voltage swing across the varactor is required if more accurate results are wished for. The circuit simulation software *Advanced Design System* provides a solution for this using a built-in component called Symbolically Defined Devices (SDD). The SDD enables the creation of equation-based nonlinear components that are a function of the port voltages, currents, and their derivates. In the figure below an alternative schematic using the SDD can be seen.

![Nonlinear SPICE model of varactor diode](https://i.imgur.com/549EoOt.png "Nonlinear SPICE model of varactor diode")

To model the nonlinear varactor using the SDD, an equation describing the capacitance as a function of reverse voltage $C(v)$ is required. This can either be obtained by fitting a curve to measured data or, if available, be found in the datasheet of the varactor. Integrating $C(v)$ with respect to $v$ gives the charge $Q(v)$ as a function of reverse voltage which can be used as input for the SDD.

The Varactor Diode

A Transmit Receive Module carries out the necessary functions in a phased array just prior to the signal being emitted through the antenna element. It encompasses both transmission and reception capabilities, which is why it’s called a Transmit Receive Module.

**A Transmit Receive Module carries out the necessary functions in a phased array just prior to the signal being emitted through the antenna element. It encompasses both transmission and reception capabilities, which is why it’s called a Transmit Receive Module.**

A Transmit Receive Module's (TRM) main components are amplifiers, filters, phase shifters, switches, and a circulator. A schematic for a typical TRM is shown in the figure. The TRM can be switched into either be in transmit or receive mode which alters the path of the signal. Some of the functionality of the TRM is used both while transmitting and receiving, namely the phase shifter and attenuator. The phase shifter is what enables the beamforming to scan the beam in different directions while the attenuator is used to control the signal amplitude. Sharing these components between both modes has several benefits. Not only does it reduce the number of active components that would increase the complexity but it also saves weight and volume. While the savings for a single TRM might be low, the total decrease in weight and volume for an array of thousands of TRMs is significant.

## Power Amplifiers in TRMs
Inside the TRM there are two main amplifiers which are shown in the figure above. In the receive path sits the low noise amplifier (LNA) that is in charge of amplifying the received signal from the antenna. This very weak signal is greatly affected by noise which is why minimizing the noise that the amplifier adds is of high importance. In addition to the LNA, a receive protection circuit is used to make sure that signals with too high power do not get into the LNA. 

In the transmit path sits the HPA. It is responsible for amplifying the signal to the desired power level and ideally provides perfectly linear amplification. However, to maximize the transmitted power a trade off between linearity and output power usually has to be made. For example, having too much output power will start to compress the output signal causing distortions and these are things an HPA designer has to take into account. Furthermore, an HPA is designed to be presented to a specific impedance, both on the input and output ports. Typically both these are set to 50 $\Omega$ by design, requiring the rest of the system to be designed around the same impedance but this varies for different applications.

In the case of an HPA with a 50 $\Omega$ output impedance that is presented to a load differing from the designed one, reflections will occur. These reflections can greatly affect the HPA resulting in reduced power efficiency and gain which can have a severe impact on its performance or potentially even break it. Making sure that the HPA is protected from reflections is therefore of high importance to achieve greater and more stable performance.

## The TRM's RF Circulator

Between the TRM and antenna, a circulator is most commonly used. A circulator is a passive component and a type of duplexer, allowing for bidirectional communication of a single path. In practice this means that signals coming from the HPA in the transmit path will be sent to the antenna and signals coming from the antenna will be sent to the receive path, thereby separating the signals without using a switch to select one path. Additionally, it also makes sure that reflections occurring at the antenna do not travel back to the HPA, degrading its performance. While the use of a circulator handles the main problems of connecting the TRM and antenna it also comes with the drawbacks of being expensive and requiring a considerable amount of volume. Another drawback is that while it typically adds less than 1 dB loss of signal power, the total loss of power can be larger when accounting for reflections between the circulator and antenna. 

## Transmit Receive Modules in Phased Antenna Arrays

The enabling technology behind Active Electronically Scanned Arrays (AESAs) comes from the ability to control hundreds or thousands of TRMs individually. Each TRM connects to an individual antenna element and is collectively controlled by a beamformer that decides how every TRM should behave at all times. This allows an AESA to create multiple directed beams which can be used for multiple tasks at once or to improve one further. 

The most common implementation of AESAs is for radar systems. Radar systems based on AESAs have several advantages compared to other technologies such as mechanically scanned arrays (MSAs) and passive electronically scanned arrays (PESAs). Both AESAs and PESAs have quicker scanning speeds than MSAs since they do not require any mechanically moving parts. 

Furthermore, the main difference between PESAs and AESAs is the number of TRMs. While a PESA has individual phase shifters, it only has a single TRM compared to AESAs which have one for each element. This gives many advantages such as allowing for the multibeam functionality of AESAs that allow radar systems to track a larger amount of targets. The individual beams can also be used to transmit multiple frequencies at once to disguise their signals which is a common task within electronic warfare and is not possible with PESAs that only can transmit a single beam. Another key feature of having individual TRMs is that even if a few TRMs fail, the functionality of the array remains which adds redundancy to the system. 

## Beamformer
Connecting to all the TRMs is the beamformer. The beamformer is what distributes and combines the signal between the TRMs in addition to steering the TRM's respective phase shift and amplitude at all times. When in transmission mode, the beamformer takes a modulated waveform and distributes it between the TRMs, and decides which phase shift should be added to each of them along with the amplitude to get the desired beam. In receive mode, the beamformer still controls the TRMs but instead of distributing a waveform, it is responsible for combining the received TRM signals into a single signal. Both in transmit and receive mode, the beamformer's phase control of the TRMs is crucial to get a beam that transmits or receives in the desired direction.

What is a Transmit Receive Module?

In this post we will run through regular white noise and iid noise and how to identify them in your residuals after applying for example a ARMA or GARCH model.

**In time series sample data there is always noise. But there exist different types of noise. In this post we will run through regular white noise and iid noise and how to identify them in your residuals after for example applying an ARMA or GARCH model. First we need to define weak stationarity before we define white noise.**

## Definition of Weak Stationarity

A time series $X$ is weakly stationary if

- $\mathbb{E}[X_t] = \mu = \text{const}$ for all $t \in \mathbb{Z}$
- $\text{Var}[X_t] = \sigma^2 = \text{const} < \infty$ for all $t \in \mathbb{Z}$
- $\text{Cov}[X_r, X_s] = \text{Cov}[X_{r+h}, X_{s+h}]$ for all $s,r,h \in \mathbb{Z}$

## Defintion of White Noise

A time series $X$ is said to be white noise if

- $X$ is weakly stationary
- $\gamma(h) = \text{Cov}[X_t, X_{t+h}] = 0$ where $h \neq 0$

## Definition of IID Noise

What is iid noise? A time series $X$ is said to be iid noise if it meets the conditions for regular white noise above but also $X_r$ and $X_s$ ($r \neq s$) is independent of each other and identically distributed for all $r,s \in \mathbb{Z}$, that is why it is called iid noise where iid stands for independent and identically distributed. 

Before going through the tricks to identifying the two different kind of noises we need to go through the definitions of the autocorrelation function (ACF) and the partial autocorrelation function (PACF).

## The Autocorrelation Function (ACF)

The autocorrelation function for the time series $X$ is defined as:

$$
    \rho_X(h) = \frac{\gamma(h)}{\gamma(0)} = \frac{\text{Cov}[X_t, X_{t+h}]}{\text{Var}[X_t]}.
$$

## The Partial Autocorrelation Function (ACF)

The autocorrelation function $\alpha(h)$ for the time series $X$ is defined by $\alpha(0) = 1$ and $\alpha(h) = \phi_{hh}$ for $h \geq 1$ where $\phi_{hh}$ is the last component of

$$
    \phi_h = ([\gamma(i-j)]_{i,j=1}^{h})^{-1} [\gamma(1), ... , \gamma(h)]'.
$$

## How to Identify White Noise

The first simple thing you could do to see if your data is just white noise is if it looks like it has no structure. The second thing is looking at the sample ACF and PACF graphs for the data. Lets look at the definition of white noise and the ACF. We know that $\gamma(h) = 0$ for $h \geq 1$ for white noise which means $\rho_X(h) = 0$ for $h \geq 1$. So all lags greater than zero should be within the red striped confidence band ($\pm 1.96/\sqrt{n}$, where $n$ is the sample size) of the ACF in the plot below. So the following sample ACF plot indicates white noise:

![Sample ACF of the sample data](https://i.imgur.com/3exyzTL.png "Sample ACF of the sample data")

The same can be said for the PACF. If $X$ is white noise then the vector $[\gamma(1), ... , \gamma(h)]'$ will just be filled with zeros and thus $\alpha(h) = 0$ for $h \geq 1$.

## How to Identify IID Noise

To identify iid noise we will have to look at the sample ACF and PACF of the data squared, $X^2$. Lets start with the ACF. We have:

$$
    \gamma_{X^2}(h) = \text{Cov}[X^2_t, X^2_{t+h}] = \mathbb{E}[(X_t^2 - \mathbb{E}[X_t^2])(X_{t+h}^2 - \mathbb{E}[X_{t+h}^2])]
$$

and now we use the property of iid noise that $X_t$ and $X_{t+h}$ is independent of each other so we get the following:

$$
    (\mathbb{E}[X_t^2 - \mathbb{E}[X_t^2]]) (\mathbb{E}[X_{t+h}^2 - \mathbb{E}[X_{t+h}^2]]) = (\mathbb{E}[X_t^2] - \mathbb{E}[\mathbb{E}[X_t^2]]) (\mathbb{E}[X_{t+h}^2] - \mathbb{E}[\mathbb{E}[X_{t+h}^2]]).
$$

Lastly we use the general property that $\mathbb{E}[\mathbb{E}[Y]] = \mathbb{E}[Y]$ so we get:

$$
    (\mathbb{E}[X_t^2] - \mathbb{E}[X_t^2]) (\mathbb{E}[X_{t+h}^2] - \mathbb{E}[X_{t+h}^2]) = 0.
$$

So the sample ACF for $X^2$ should be in the confidence band for all lags greater than zero. For the PACF we can see that the vector $[\gamma_{X^2}(1), ... , \gamma_{X^2}(h)]'$ would be zero for all $h$. Therefore the PACF should also be within the confidence band for lags greater than zero if $X$ were to be iid noise which we can see in the following plot:

![Sample PACF of the sample data squared](https://i.imgur.com/4ywfylN.png "Sample ACF of the sample data squared")

So to conclude, for iid noise $X$ and $X^2$'s sample ACF and PACF should all have to be within the confidence band for lags greater than zero.

## Ljung-Box Test for IDD Noise

Lastly one can do a Ljung-Box test to determine whether the data is independently distributed. The test is defined by the following hypothesis:

$$
        \begin{aligned}
            & H_0 : X \sim \text{IID}(\mu, \sigma^2) \\
            & H_1 : X \not\sim \text{IID}(\mu, \sigma^2) 
        \end{aligned}
$$

The Ljung-Box test has a test statistic which is:

$$
    \lambda = n (n + 2) \sum_{i=1}^h \frac{\hat{\rho}(i)^2}{n-i},
$$

where $\hat{\rho}(i)$ is the sample ACF at lag $i$, $n$ is the sample size and $h$ is the number of lags being tested. The null-hypothesis, $H_0$, is rejected at the $\alpha$-level, $\alpha \in (0,1)$, if $\lambda > \chi_{1-\alpha, h}^2$, where $\chi_{1-\alpha, h}^2$ denotes the $1 - \alpha$-quantile of the $\chi^2$-distribution with $h$ degrees of freedom.

Time Series White Noise

The continuity equation is an important tool that one can use whether it is in a fluid dynamics problem or an electromagnetic problem. The continuity equation for charges relates the charge density to the current density in a volume.

**The continuity equation is an important tool that one can use whether it is in a fluid dynamics problem or an electromagnetic problem. The continuity equation for charges relates the charge density to the current density in a volume.**

## Derivation of Continuity Equation for Charges

Lets derive the continuity equation for electric charges. If you are just searching for it, this is what it looks like:

$$
    \frac{\partial \rho}{\partial t} = - \nabla \cdot \vec{j},
$$

where $\rho$ is the charge density and $\vec{j}$ is the current density (current divided by an area). Lets get a few things straight first. When we say the charge is conserved we are talking about the charge in an arbitrary volume. Additional charge can not be created from some sort of source inside the volume. The only way to change the charge is a flow of new charges from outside the volume, that would be our $\vec{j}$. So lets observe an infinitesimal surface element $d \vec{S}$ on the surface of our volume. The surface element is vector that is normal to the surface. So the infinitesimal charge $d q$ that passes through the surface element can be written as:

$$
    dq = \vec{j} \cdot d \vec{S} dt.
$$

We'll use this expression later. The total charge $q$ for the whole volume $V$ can be written as:

$$
    q = \int_V \rho dV,
$$

where earlier mentioned $\rho$ is the charge density for the volume. We can now use the conservation of charge and look at the change of the charge through time and flux of charge through the volume's surface $\partial V$. So we have two expressions, firstly:

$$
    \frac{dq}{dt} = \int_V \frac{\partial \rho}{\partial t} dV,
$$

which is just the derivative of the total charge for the volume. Secondly, we have:

$$
    \frac{dq}{dt} = - \int_{\partial V} \vec{j} \cdot d \vec{S},
$$

which uses the expression for the charge through a surface element. The reason for the minus sign is because we are calculating the inflow of charges but the surface element vector is pointing outwards from the volume and we have to counteract that. We are almost at our sought after expression but we'll have to use one more trick, the divergence theorem aka Guass's theorem which states:

$$
    \int_V \nabla \cdot \vec{F} dV = \int_{\partial V} \vec{F} \cdot d \vec{S},
$$

which fits perfectly with our charge flux expression:

$$
     - \int_{\partial V} \vec{j} \cdot d \vec{S} = - \int_V \nabla \cdot \vec{j} dV.
$$

So now we have:

$$
    \frac{dq}{dt} = \int_V \frac{\partial \rho}{\partial t} dV = - \int_V \nabla \cdot \vec{j} dV.
$$

As you can see we have two expression with a volume integral for the volume $V$. So we can remove the integral and we finally get our continuity equation for electric charges:

$$
    \frac{\partial \rho}{\partial t} = - \nabla \cdot \vec{j}.
$$

Now lets see how the continuity equation looks in different coordinate systems.

## Continuity Equation in Cartesian Coordinates

The continuity equation for charges in Cartesian coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{\partial j_x}{\partial x} + \frac{\partial j_y}{\partial y} + \frac{\partial j_z}{\partial z} = 0.
$$

## Continuity Equation in Cylindrical Coordinates

The continuity equation for charges in cylindrical coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r}(r j_r) + \frac{1}{r} \frac{\partial}{\partial \theta}(j_\theta) + \frac{\partial}{\partial z}(j_z) = 0.
$$

## Continuity Equation in Spherical Coordinates

The continuity equation for charges in spherical coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 j_r) + \frac{1}{r \sin{\theta}} \frac{\partial}{\partial \theta}(j_\theta \sin{\theta}) + \frac{1}{r \sin{\theta}} \frac{\partial}{\partial \varphi}(j_\varphi) = 0.
$$

Continuity Equation in Electromagnetism

The moment of inertia is dependent on where the axis of rotation is located. Lets calculate the moment of inertia for the different rotational axes for a rod and a disk.

**If you are just searching for the moment of inertia of a rod or a disk I will just give them to you. For a rod $I = \frac{L^2}{12}$ if the rotation axis is located at the middle the rod and $I = M \frac{L^2}{3}$ if the rotational axis is located at the end of the rod. For a homogeneous disk/cylinder $I = M \frac{R^2}{2}$ and for a disk/cylinder with just a shell $I = MR^2$. Now if you are interested in how we get these values feel free to read further.**

We know that the moment of a point of inertia of a point mass is

$$
    M R^2,
$$

where $M$ is the mass and $R$ is the distance between the mass and the rotational axis. We can use this to calculate the moment of inertia for objects that is not point masses. We do this by integrating over all infinitesmal mass elements $dm$ of the object. So our moment of inertia is now

$$
    I = \int_0^M r^2 dm.
$$

## The Moment of Inertia of a Rod

The moment of inertia is dependent on the location of the rotational axis, therefore the moment of inertia will vary for a rod depending on the axis of rotation. We have our general formula for the moment of inertia

$$
    I = \int_0^M r^2 dm.
$$

We now have to identify what $dm$ is. Since the rod can be looked at in one dimension

$$
    dm = \frac{M}{L} dr
$$

where $L$ is the length of the rod. $\frac{M}{L}$ can be looked at as an one dimensional density for the rod. Now lets calculate the moment of inertia if the rod is rotating at one of its ends. The moment of inertia then becomes

$$
    I = \int_0^L r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_0^L = \frac{M}{L} \frac{L^3}{3} = M \frac{L^2}{3}.
$$

If the axis of rotation is in the middle of the rod the moment of inertia then becomes

$$
    I = \int_{-L/2}^{L/2} r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_{-L/2}^{L/2} = \frac{M}{3L} \left[ \frac{L^3}{8} - \left(\frac{-L^3}{8} \right) \right] = M \frac{L^2}{12}.
$$

Lets create a more general formula. Lets use the variable $x$ for the point of rotation. Our moment of inertia integral then becomes

$$
    I = \int_{-x}^{L-x} r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_{-x}^{L-x} = \frac{M}{3L} \left[ (L-x)^3 + x^3 \right].
$$

## The Moment of Inertia of a Disk

![Image of the moment of inertia for disks](https://i.imgur.com/lyfEGNK.png "Moment of inertia for a homogeneous disk and a hollow disk")

Lets move on to calculate the moment of inertia for a disk (this also applies to a cylinder which is kind of an elongated disk). The axis of rotation will be located in the center of the disk but we will get back to calculating the moment of inertia of a displaced rotational axis. Like usual we start of with the classic integral for the moment of inertia

$$
    I = \int_0^M r^2 dm.
$$

Defining $dm$ is a little bit more tricky than for the rod. In this problem we have to go two dimensional. Lets start with defining our two dimensional density which is

$$
    \frac{M}{\pi R^2}.
$$

The denominator is just the area of the disk. Our area element is 

$$
    r dr d\theta
$$

I won't go into to much detail here, you'll have to do some repetition on polar coordinates if don't understand how we get this area element. So our moment of inertia becomes

$$
    I = \frac{M}{\pi R^2} \int_0^{2 \pi} \int_0^R r^3 dr d\theta = \frac{2M}{R^2} \int_0^R r^3 dr = \frac{2M}{R^2} \Bigg[ \frac{r^4}{4} \Bigg]_0^R = \frac{2M}{R^2} \frac{R^4}{4} = M \frac{R^2}{2}.
$$

Now, what happens if the axis of rotation was shifted (but still parallel to the axis through the center of mass)? Well luckily for us we don't have to do a super complicated integral, we can just use the parallel axis thereom which says

$$
    I = I_\text{CM} + Md^2,
$$

where $I_\text{CM}$ is the moment of inertia through the center of mass, in this case $M \frac{R^2}{2}$. $d$ is the distance between the axis of rotation and the rotational axis throught the center of mass. You can test this formula by calculating the moment of inertia at an end of the rod. In that case we have $I_\text{CM} = M \frac{L^2}{12}$ and $d=\frac{L}{2}$.

So what if the disk/cylinder is hollow? Then we simply change the limits of integration and the density. If the disk is hollow up to $r=a$ and then has mass up to $r=b$ the 2D density changes to

$$
    \frac{M}{\pi (b^2-a^2)}.
$$

So we get the following moment of inertia

$$
    I = \frac{2M}{b^2-a^2} \int_a^b r^3 dr = \frac{2M}{b^2-a^2} \Bigg[ \frac{r^4}{4} \Bigg]_a^b = \frac{2M}{b^2-a^2} \frac{b^4-a^4}{4} = M \frac{b^2+a^2}{2}.
$$

Here we can see that if the disk/cylinder only consisted of a shell ($a \rightarrow b$) we get

$$
    I = M b^2.
$$

So now you know how to calculate the moment of inertia of a rod and a disk.

The Moment of Inertia of a Rod and a Disk

What makes the suspension wire of the Golden Gate Bridge curve like a parabola? Lets get into the mechanics of the shape of the suspension bridge's cables and what determines them to curve like a parabola or a catenary.

**What makes the suspension wire of the Golden Gate Bridge curve like a parabola? Lets get into the mechanics of the shape of the suspension bridge's cables and what determines them to curve like a parabola or a catenary.**

## The Assumptions about the Problem

Some assumptions about the problem have to be made first. The main wire that the bridge rests on is assumed to be inextensible, homogeneous and infinitely flexible. The bridge is also assumed to be homogeneous. The main wire is connected at each end of the bridge. The wires that ties together the bridge and the main wire is assumed to be in an infinite amount so the bridge's load function is constant. The origin of coordinate system is located at the bottom of the main wire.

## Equilibrium of an Infinitesimal Wire Element

Lets make an free body diagram of an infinitesmal wire element of the main wire:

![Image of free body diagram of an infinitesmal wire element](https://i.imgur.com/vGGtqIv.png "Free body diagram of an infinitesmal wire element")

Here $\vec{T}$ is the tension of the cable, $d\vec{T}$ the infinitesmal change in the tension, $ds$ the infinitesmal length of the wire element and $dx$ the length in the $x$-direction. $q_c$ and $q_b$ is respectively the cable's and bridge's linear mass density. The equilibrium of the forces is

$$
    \vec{T} + d\vec{T} - \vec{T} - (q_c ds + q_b dx)g \hat{y} = 0.
$$

With some cancellation we get

$$
    d\vec{T} = (q_c ds + q_b dx)g \hat{y}.
$$

Lastly, the Pythagorean theorem gives us 

$$
    ds = \sqrt{(dx)^2 + (dy)^2} = \sqrt{1 + (y')^2} dx.
$$

So we get

$$
    d\vec{T} = (q_c \sqrt{1 + (y')^2} dx + q_b dx)g \hat{y}.
$$

## The ODE which Dictates the Shape of The Wire

The tension at the bottom of the main wire has only a $x$-component which we'll call $T_0$. So $\vec{T}(0) = T_0 \hat{x}$. If we do some integration of the equilibrium equation from earlier we get

$$
    \vec{T} = T_0 \hat{x} + \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right)g \hat{y},
$$

here we are just using $\tilde{x}$ in the integral so the upper limit, $x$, of the integral is correct. In the $x$-direction we now have

$$
    T \cos{\theta} = T \frac{dx}{\sqrt{(dx)^2 + (dy)^2}} = T_0,
$$

where $\theta$ is the angle between $\vec{T}$ and the $x$-axis. So

$$
    T = T_0 \frac{\sqrt{(dx)^2 + (dy)^2}}{dx}.
$$

In the $y$-direction we have

$$
    T \sin{\theta} = T \frac{dy}{\sqrt{(dx)^2 + (dy)^2}} = T_0 \frac{dy}{dx} = \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right)g.
$$

We now have an expression for $y'$:

$$
    y' = \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right) \frac{g}{T_0}
$$

Lets do an additional derivation to get rid of the integral:

$$
    y'' = \frac{q_b g}{T_0} + \frac{q_c g \sqrt{1+(y')^2}}{T_0}.
$$

This is the differential equation that dictates the curve of suspension bridge.

## Is the Cable Shaped as a Parabola or Catenary?

So what shape is the cable? To answer this lets look at two extremes, neglecting either the cable's mass or the bridge's mass. If we neglect the cable's load function i.e. its mass we get the differential equation

$$
    y'' = \frac{q_b g}{T_0}.
$$

With some simple integration we get

$$
    \frac{q_b g}{2 T_0} x^2 + bx + c,
$$

where $b$ and $c$ are constants. The origin is loacted at the bottom of the cable so $c=0$ and the cable is symmetric so $b=0$. Here we have our parabola. Neglecting the cable's mass is a pretty reasonable approximation since the bridge's mass is likely much larger than the cable's mass. Can you guess what we get if we do the opposite? If we neglect the bridge's mass we get the following differential equation:

$$
    y'' = \frac{q_c g \sqrt{1+(y')^2}}{T_0}.
$$

Since we already have pretty good hunch that the solution is in the shape of a catenary, lets check if $y = a \cosh{\frac{x}{a}}$ works. Then the left side is

$$
    y'' = \frac{\cosh{\frac{x}{a}}}{a}
$$

and the right side is

$$
    \frac{q_c g \sqrt{1+(y')^2}}{T_0} = \frac{q_c g \sqrt{1+\sinh^2{\frac{x}{a}}}}{T_0} =  \frac{q_c g \cosh{\frac{x}{a}}}{T_0}.
$$

For this to work $a$ needs to be

$$
    a = \frac{T_0}{q_c g}.
$$

That is the solution if the bridge's mass can be neglected. I don't know what sort of bridge that is... But hey now you know. I'll leave the harder task of solving the ODE with both masses present for you.

Is the Suspension Bridge Shaped like a Parabola or a Catenary?

When studying transient processes, like when a switch is changed, in electrical circuits we have to make us of the Laplace transform. Here is a quick rundown on how to Laplace transform the resistor, capacitor and inductor.

**When studying transient processes, like when a switch is changed, in electrical circuits we have to make us of the Laplace transform. Here is a quick rundown on how to Laplace transform the resistor, capacitor and inductor.**

First we need to establish what happens when we Laplace transform the current and voltage. We can see the voltage and current as "unknown" functions dependent on time, therefore their Laplace transforms become:

$$
        \begin{aligned}
            & i(t) \Rightarrow I(s) \\
            & u(t) \Rightarrow U(s)
        \end{aligned}
        .
$$

## Laplace Transform of the Resistor

So lets start with the easy component: the resistor. We know from Ohm's law that:

$$
    u(t) = R i(t).
$$

If we Laplace transform this we get:

$$
    U(s) = R I(s),
$$

since the resistor is time independent there is not much happening when moving from the t-domain to the s-domain. The resistor remains the same in both domains.

## Laplace Transform of a Derivative

The inductor and capacitor is time dependent, so they get a little bit more trickier. First we need to see what happens when we apply the Laplace transform to a derivative of a function $f(t)$. The Laplace transform of an arbitrary function $f(t)$ is given by:

$$
    F(s) = \int_{0-}^\infty e^{-st} f(t) dt.
$$

The $0-$ can be looked at more or less as what happened exactly before significant change happened to the function $f(t)$ (we'll get back to this). So lets preform the Laplace transform on $f'(t)$:

$$
    \int_{0-}^\infty e^{-st} f'(t) dt.
$$

With a little partial integration we get:

$$
    \bigg[e^{-st} f(t) \bigg]_{0-}^\infty - \int_{0-}^\infty -se^{-st} f(t) dt = - f(0-) + sF(s) .
$$

## Laplace transform of the Inductor

With that out of the way lets preform the Laplace transform on the inductor first. We know that the potential over an inductor is given by:

$$
    u(t) = L \frac{d i(t)}{dt}.
$$

So with our knowledge of Laplace transforming a derivative we get the following:

$$
    U(s) = L (s I(s) - i(0-)).
$$

The inductor is said to be energized if there is an initial current $i(0-)$. If there is no initial current we get a simpler relation of:

$$
    U(s) = sL I(s).
$$

We can also rewrite the relation to get the current:

$$
    I(s) = \frac{U(s)}{sL} + \frac{i(0-)}{s}.
$$

We can make circuit diagrams with these relations:

![Image of Laplace transform of the inductor](https://i.imgur.com/IkgvQUe.png "Laplace transform of the inductor")

Note the Thevenin/Norton equivalents under "Initial current".

## Laplace transform of the Capacitor

Lets move on to the capacitor. The procedure of Laplace transforming the capacitor is the same as the inductor. The current through a capacitor is given by:

$$
    i(t) = C \frac{du(t)}{dt}.
$$

So the Laplace transform of the current is:

$$
    I(s) = C(sU(s) - u(0-)).
$$

We can see here that if there is no initial voltage over the capacitor we get the simpler relation:

$$
    I(s) = sC U(s).
$$

Lastly the the Laplace transformed voltage is:

$$
    U(s) = \frac{I(s)}{sC} + \frac{u(0-)}{s}.
$$

The circuit diagrams then become:

![Image of Laplace transform of the capacitor](https://i.imgur.com/za9EeTS.png "Laplace transform of the capacitor")

And that is it! You now know how to Laplace transform the resistor, inductor and capacitor!

Laplace Transform of DC Circuits

The Navier Stokes equation is famously known for being unsolvable if there is any turbulence in a fluid but there are certain scenarios where laminar flow exist and a solution for the velocity profile can be derived. One of these scenarios is flow past two plates either caused by one plate moving, Couette flow, or a pressure gradient.

**The Navier Stokes equation is famously known for being unsolvable if there is any turbulence in a fluid but there are certain scenarios where laminar flow exist and a solution for the velocity profile can be derived. One of these scenarios is flow past two plates either caused by one plate moving, Couette flow, or a pressure gradient.**

Before solving these two problems we need to make some assumptions. First we assume that our fluid is incompressible and Newtonian so we can use the Navier Stokes equation and the mass continuity equation. Secondly, the flow has to be stationary (so all time dependent terms in the Navier Stokes equation is zero) and in one direction, the $x$-direction, so the velocity vector only has one component ($v=w=0$). The plates also has to be considered infinite. Lastly, the flow has to be laminar which has already been mentioned and the impact from gravity has to be considered negligible ($g = 0$).

## Couette Flow - One Plate Moves

So lets start with the case where one of the plates is moving. From our assumptions that fluid is Newtonian, stationary and goes in the $x$-direction we get the following mass continuity equation:

$$
    \frac{\partial u}{\partial x} + \cancel{\frac{\partial v}{\partial y}} + \cancel{\frac{\partial w}{\partial z}} = \frac{\partial u}{\partial x} = 0,
$$

where $u$, $v$ and $w$ is the velocity in the $x$-, $y$- and $z$-direction. It is then safe to assume that the velocity profile $u$ is only dependent on $y$ because the plates' areas are infinite. Now lets look at the Navier Stoke equation in all directions and see which terms we can cancel out because of our assumptions and the mass contuinity equation:

$$
    \left\{
        \begin{aligned}
            & \rho \left( \cancel{\frac{\partial u}{\partial t}} + \cancel{u \frac{\partial u}{\partial x}} + \cancel{v \frac{\partial u}{\partial y}} + \cancel{w \frac{\partial u}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial x}} + \mu \left(\cancel{\frac{\partial^2 u}{\partial x^2}} + \frac{\partial^2 u}{\partial y^2} + \cancel{\frac{\partial^2 u}{\partial z^2}} \right) + \cancel{\rho g_x} \\
            & \rho \left( \cancel{\frac{\partial v}{\partial t}} + \cancel{u \frac{\partial v}{ \partial x}} + \cancel{v \frac{\partial v}{\partial y}} + \cancel{w \frac{\partial v}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial y}} + \mu \left(\cancel{\frac{\partial^2 v}{\partial x^2}} + \cancel{\frac{\partial^2 v}{\partial y^2}} + \cancel{\frac{\partial^2 v}{\partial z^2}} \right) + \cancel{\rho g_y} \\
            & \rho \left( \cancel{\frac{\partial w}{\partial t}} + \cancel{u \frac{\partial w}{ \partial x}} + \cancel{v \frac{\partial w}{\partial y}} + \cancel{w \frac{\partial w}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial z}} + \mu \left(\cancel{\frac{\partial^2 w}{\partial x^2}} + \cancel{\frac{\partial^2 w}{\partial y^2}} + \cancel{\frac{\partial^2 w}{\partial z^2}} \right) + \cancel{\rho g_z}
        \end{aligned}
    \right. .
$$

As you can see, a bunch of terms cancel out and all that remains is just one single expression:

$$
    \mu \frac{\partial^2 u}{\partial y^2} = 0 \Rightarrow \frac{\partial^2 u}{\partial y^2} = 0.
$$

So if we integrate this expression twice we get:

$$
    u(y) = C_1 y + C_2,
$$

where $C_1$ and $C_2$ are constants. To determine these constants we have to use the no-slip boundary conditions of the plates. If we assume that the origin is located exactly in between the two plates we know that:

$$
    \left\{
        \begin{aligned}
            & u(+h) = u_p \\
            & u(-h) = 0
        \end{aligned}
    \right. ,
$$

where $h$ is the shortest distance from the origin to a plate and $u_p$ is the speed of the moving plate. This gives the constants:

$$
    \left\{
        \begin{aligned}
            & C_1 = \frac{u_p}{2 h} \\
            & C_2 = \frac{u_p}{2}
        \end{aligned}
    \right. .
$$

So our velocity profile becomes:

$$
    u(y) = \frac{u_p}{2} \left(1 + \frac{y}{h} \right).
$$

This solution is called the Couette flow which is named after the french physicist Maurice Couette. 

## Laminar Flow Past Two Plates Caused by Pressure Gradient

So what solution do we get if the two plates are stationary and we instead have a constant pressure gradient? The assumptions mentioned earlier are still the same so we get the same mass continuity equation. So lets see what cancels out in our Navier Stokes equations:

$$
    \left\{
        \begin{aligned}
            & \rho \left( \cancel{\frac{\partial u}{\partial t}} + \cancel{u \frac{\partial u}{\partial x}} + \cancel{v \frac{\partial u}{\partial y}} + \cancel{w \frac{\partial u}{\partial z}} \right) = - \frac{\partial p}{\partial x} + \mu \left(\cancel{\frac{\partial^2 u}{\partial x^2}} + \frac{\partial^2 u}{\partial y^2} + \cancel{\frac{\partial^2 u}{\partial z^2}} \right) + \cancel{\rho g_x} \\
            & \rho \left( \cancel{\frac{\partial v}{\partial t}} + \cancel{u \frac{\partial v}{ \partial x}} + \cancel{v \frac{\partial v}{\partial y}} + \cancel{w \frac{\partial v}{\partial z}} \right) = - \frac{\partial p}{\partial y} + \mu \left(\cancel{\frac{\partial^2 v}{\partial x^2}} + \cancel{\frac{\partial^2 v}{\partial y^2}} + \cancel{\frac{\partial^2 v}{\partial z^2}} \right) + \cancel{\rho g_y} \\
            & \rho \left( \cancel{\frac{\partial w}{\partial t}} + \cancel{u \frac{\partial w}{ \partial x}} + \cancel{v \frac{\partial w}{\partial y}} + \cancel{w \frac{\partial w}{\partial z}} \right) = - \frac{\partial p}{\partial z} + \mu \left(\cancel{\frac{\partial^2 w}{\partial x^2}} + \cancel{\frac{\partial^2 w}{\partial y^2}} + \cancel{\frac{\partial^2 w}{\partial z^2}} \right) + \cancel{\rho g_z}
        \end{aligned}
    \right. .
$$

We can see that the pressure gradient is zero in the $y$- and $z$-direction. So we get the following differential equation:

$$
    \mu \frac{\partial^2 u}{\partial y^2} = \frac{\partial p}{\partial x}.
$$

With some integration we get:

$$
    u(y) = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{y^2}{2} + C_1 y + C_2,
$$

where $C_1$ and $C_2$ are constants. The velocity at both plates are now zero because of the no-slip condition, so:

$$
    \left\{
        \begin{aligned}
            & u(+h) = 0 \\
            & u(-h) = 0
        \end{aligned}
    \right. .
$$

This gives the constants:

$$
    \left\{
        \begin{aligned}
            & C_1 = 0 \\
            & C_2 = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{h^2}{2}
        \end{aligned}
    \right. .
$$

So our velocity profile becomes:

$$
    u(y) = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{h^2}{2} \left(1 + \left(\frac{y}{h} \right)^2 \right).
$$

The velocity profile is a parabola with its maximum at the middle of the plates. 

![Image of pressure gradient flow](https://i.imgur.com/xIiH6ge.png "Pressure Gradient Flow")

So that is it, that is how you calculate fluid flow past two plates!

Couette Flow - Laminar Flow Past Two Plates

Why is it that a falling chimney breaks at a third of its length? Lets find out why!

**The falling chimney problem is the perfect mix of mechanical dynamics and beam theory. By just using the fundamentals of those two fields we can accurately predict where the chimney cracks when falling, under some assumptions...**

So lets get the assumptions of the falling chimney problem out of the way first. The first assumption is that there is no air resistance. If air resistance had been present the top of the chimney would experience a greater stress from the air than the lower part which makes our calculations much more complex and this is not an article about fluid mechanics. The second assumption is that the chimney is homogeneous, meaning that it has the same diameter across its length and the same material. If it had been inhomogeneous certain parts would handle stress better than others.

## Angle Acceleration of the Chimney

Now when that is out of the way, lets make some calculations! The first thing we'll have to determine is the angle acceleration $\ddot{\theta}$. To do that we use a beam's, in this case a chimney's, moment of inertia when attached at one of its ends:

$$
    I = \frac{1}{3} m L^2.
$$

Here $m$ is the chimneys mass and $L$ is its length. The moment from the gravitation acting on the chimney must equal that of the angle acceleration:

$$
    mg \frac{L}{2} \sin{\theta} = I \ddot{\theta} = \frac{1}{3}mL^2 \ddot{\theta}.
$$

So now we have our angle acceleration $\ddot{\theta}$ and apparently it is dependent on $\theta$:

$$
    \ddot{\theta} = \frac{3g}{2L} \sin{\theta}.
$$

Is this plausible? Well $\ddot{\theta}$ increases with $g$ and decreases with $L$ which makes sense. At $\theta=0$ the angle acceleration is zero which also makes sense since at that angle the chimney is "standing up".

## The Chimney's Load Function

If the chimney was locked and stationary at an angle $\theta$ the load function would be:

$$
    w(d) = -\frac{mg}{L} \sin{\theta}.
$$

If we go to the other extreme and say that the chimney was in free fall it would lead to the load function being zero. So falling negates the load function and in our case the chimney is rotating. So $\ddot{\theta}$ has a similar effect to free fall. So the load function in our case will be the following:

$$
    w(d)=\frac{dm}{L} \ddot{\theta} - \frac{mg}{L} \sin{\theta} = \frac{3dmg}{2L^2} \sin{\theta} - \frac{mg}{L} \sin{\theta}.
$$

We can see here that the load function is zero at $d = \frac{2L}{3}$. This makes sense because $\ddot{\theta}$ is equal to $g \sin{\theta}$ at $d = \frac{2L}{3}$. In other words, this point is in a "virtual" free fall.

## The Chimney's Shear Force Function

The relation between the load function and the shear force function is:

$$
    -w(d)=V'(d).
$$

If we integrate that we get:

$$
    V(d) = -\frac{3d^2mg}{4L^2} \sin{\theta}+ \frac{dmg}{L} \sin{\theta}+ C,
$$

where $C$ is a constant. We can determine this constant with the boundary condition $V(L) = 0$. In other words: the shear force is zero at the end of the chimney which is logical. This means:

$$
    C = -\frac{mg}{4}.
$$

If we take the absolute value of the shear force function we can check where the force is the greatest:

$$
    |V(d)|=\left|mg\, \sin{\theta} \left(-\frac{3d^2}{4L^2} + \frac{d}{L} - \frac{1}{4}\right) \right|.
$$

Here we can see that the shear force is the greatest when $d \to 0$ on the interval $0< d < L$, so if the shear force was the cause of the crack in the chimney it would be at the base. But that is not the only place where the chimney can break... 

## The Chimney's Bending Moment Function

The final step is looking at the bending moment function which has the following relation with the shear force function:

$$
    M'(d)=V(d).
$$

So with a bit of integration we get:

$$
    M(d)=mg\, \sin{\theta} \left(-\frac{d^3}{4L^2} + \frac{d^2}{2L} - \frac{d}{4}\right) + K,
$$

where $K$ is a constant. We'll once again use the boundary condition at the end of the chimney, which in this case is also zero. This leads to $K=0$. So where is $|M(d)|$ the greatest value? At $d=0$ the bending moment is zero so it is not there. Well lets check where its derivative, the shear force function, is zero to get the maximums:

$$
    M'(d)=V(d)=0 \Longrightarrow -\frac{3d^2}{4L^2} + \frac{d}{L} - \frac{1}{4} = 0.
$$

From this we get the roots $d=L$ and $d=\frac{L}{3}$. $d=L$ does not make much sense since its the end of the chimney, you cant real chop of something that isn't there... So we get $d=\frac{L}{3}$ and that is why a falling chimney breaks at a third of its length from its base.

<div style="position: relative;
    padding-bottom: 56.25%; /* 16:9 */
    height: 0;">
<iframe style="position: absolute;
    top: 0;
    left: 0;
    width: 100%;
    height: 100%;" width="560" height="315" src="https://www.youtube.com/embed/efjegZ6aWKI?start=149" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

</div>

The video above shows a demolition of a smokestack, as you can see it breaks close to a third of its length when it falls.

The Falling Chimney

The chemical potential has different characteristics depending on the temperature. Lets calculate and understand the approximations of the chemical potential in an electron doped semiconductor.

**The chemical potential has different characteristics depending on the temperature. Lets calculate and understand the approximations of the chemical potential in an electron doped semiconductor.**

## The Chemical Potential at Low Temperatures

Lets start with temperatures close to absolute zero. You have probably heard or read that the chemical potential must lie in the center between the donor level and the bottom of the conduction band because of charge neutrality. Lets do some math to find out why that is.

Lets use the variable $E_g$ for the bandgap energy, $E_d$ is the donor energy so the the donor level is situated at $E_g-E_d$. $\mu$ is our desired chemical potential. $n$ is the electron density in the conduction band, $n_d$ is the donor density and $n_d^+$ is the density of the ionized donors. From the Fermi-Dirac distribution we know that the probability of an electron at the energy $E_g-E_d$ and the temperature $T$ is:

$$
    f(E_g - E_d) = \frac{1}{e^{(E_g-E_d-\mu)/(k_B T)}+1}.
$$

So the donor electron density at the donor level is:

$$
    \frac{n_d}{e^{(E_g-E_d-\mu)/(k_B T)}+1}.
$$

Every donor electron not at the donor level is at the conduction band so:

$$
    n_d^+ = n_d (1 - \frac{1}{e^{(E_g-E_d-\mu)/(k_B T)}+1}) = \frac{n_d}{e^{-(E_g-E_d-\mu)/(k_B T)}+1}.
$$

Since we are at low temperatures $T \ll E_d/k_B$. This means $-(E_g-E_d-\mu)/(k_B T) \gg 1$ so we can make the following approximation (the 1 in the denominator is neglected):

$$
    n_d^+ \approx n_d e^{(E_g-E_d-\mu)/(k_B T)}.
$$

Having low temperatures also means that the electrons at the valence band don\'t have enough energy to make an excitation to the conduction band meaning $n \approx n_d^+$. But the electron density $n$ can also be expressed as:

$$
    n = \frac{1}{\sqrt{2}} (\frac{m_e^* k_B T}{\pi \hbar^2})^{3/2} e^{-(E_g-\mu)/(k_B T)} = n_0 e^{-(E_g-\mu)/(k_B T)}
$$

here is $m_e^*$ the effective mass of the electrons in the conduction band. Lets put the different expressions for $n$ at low temperature equal to each other to get the chemical potential:

$$
    n_d e^{(E_g-E_d-\mu)/(k_B T)} = n_0 e^{-(E_g-\mu)/(k_B T)}
$$

So now we get the chemical potential:

$$
    \mu = E_g - \frac{1}{2} E_d - \frac{1}{2} k_B T \ln{(\frac{n_0}{n_d})}.
$$

So when the temperature approaches zero, $T \to 0$, the chemical potential becomes:

$$
    \mu = E_g - \frac{1}{2}E_d,
$$

since polynomial functions dominates over logarithmic functions ($n_0$ contains $T^{3/2}$) when it comes to limits.

## The Chemical Potential at Medium Temperatures

So that is how the chemical potential behaves at low temperatures, how about "medium" temperatures like room temperature? Well here we have more thermal energy to work with so pretty much every donor atom is ionized ($n_d^+ \approx n_d$) but the thermal energy is still not big enough to excite a significant amount of electrons in the valence band so $n \approx n_d$. This makes the math easier, we now have:

$$
    n = n_0e^{-(E_g-\mu)/(k_B T)} = n_d,
$$

so:

$$
    \mu = E_g + k_B T \ln{(\frac{n_d}{n_0})},
$$

for "medium" temperatures.

## The Chemical Potential at High Temperatures

For high temperatures, like 1200 K, electron holes has manifested at the valence band and now plays the bigger role in the electron density in the conduction band. So now $n \approx p$, where $p$ is the electron hole density at the valence band. The semiconductor is said to be in an intrinsic state, meaning that it can be seen as undoped. So now we have the following relation:

$$
    \frac{1}{\sqrt{2}} (\frac{m_e^* k_B T}{\pi \hbar^2})^{3/2} e^{-(E_g-\mu)/(k_B T)} = \frac{1}{\sqrt{2}} (\frac{m_h^* k_B T}{\pi \hbar^2})^{3/2} e^{-\mu/(k_B T)},
$$

where $m_h^*$ is the effective mass of the electron holes. This leads to the chemical potential of:

$$
    \mu = \frac{1}{2} E_g + \frac{3}{4} k_B T \ln{(\frac{m_h^*}{m_e^*})}.
$$

The chemical potential is located midway between the conduction band and the valence band if $m_e^* = m_h^*$.

So there we have it! We have now calculated the chemical potential for different temperature regions for an n-type semiconductor, although we did it with some significant approximations but without those this would get a lot harder if not impossible.

The Chemical Potential of an n-type Semiconductor

Calculating the energy of the stationary state of a two electron system, like a helium atom or lithium ion, gets quiet a bit more tricky than for a single electron system. In fact, calculating the Schrödinger equation is impossible. Luckily we have methods like perturbation theory and the variational method to estimate the energy.

**Calculating the energy of the stationary state of a two electron system, like a helium atom or lithium ion, gets quiet a bit more tricky than for a single electron system. In fact, calculating the Schrödinger equation is impossible. Luckily we have methods like perturbation theory and the variational method to estimate the energy.**

First we need to establish the Hamiltonian for the system which is:
$$
    \hat{H} = -\frac{\hbar}{2 m_e} \nabla^2_1 - \frac{\hbar}{2 m_e} \nabla^2_2 - \frac{Ze^2}{4 \pi \epsilon_0 r_1 } - \frac{Ze^2}{4 \pi \epsilon_0 r_2 } + \frac{e^2}{4 \pi \epsilon_0 | \vec{r}_1 - \vec{r}_2 | }.
$$
$\frac{\hbar}{2 m_e} \nabla^2_1$ is the kinetic energy for the first electron, same goes for the second. $\frac{Ze^2}{4 \pi \epsilon_0 r_1 }$ is Coulomb's law between the atom core's charge $Ze$ and the first electron's charge $e$. Lastly we have $\frac{e^2}{4 \pi \epsilon_0 | \vec{r}_1 - \vec{r}_2 | }$ which is Coulomb's law between the two electrons, this term is positive because the electrons repel each other.

## Two Independent Electrons Estimation

The first estimation we'll do is having the two electrons not react with each other which means that the Coulomb's law term between them is zero. Then the combined wave function of the system is the following:

$$
    \psi_+ (\vec{r}_1, \vec{r}_2) = \psi_{100}(r_1) \psi_{100}(r_2),
$$

where

$$
    \psi_{100}(r) = \frac{Z^{3/2}}{\sqrt{\pi} a_0^{3/2}} e^{Zr/a_0}
$$

is the wave function of a single electron system, $a_0$ is the Bohr radius. So our wave function of the two electron system is

$$
    \psi_+ (r_1, r_2) = \frac{Z^{3}}{\sqrt{\pi} a_0^{3}} e^{Z(r_1+r_2)/a_0},
$$

which leads the stationary state energy of:

$$
    E = 2Z^2 \cdot (-13.6 \text{ eV}).
$$

## Perturbation Theory Estimation

But that is not so satisfactory nor a good estimate, which leads to our second estimation method: perturbation theory. Lets use the first order correction to estimate the impact that Coulomb's law between the two electrons has on the energy:

$$
    E^{(1)} = \langle \psi_+ | \frac{e^2}{4 \pi \epsilon_0} \frac{1}{|\vec{r}_1 - \vec{r}_2|} | \psi_+ \rangle = \frac{5}{4} Z \cdot 13.6 \text{ eV}.
$$

So now our estimated stationary state energy is:

$$
    E = 2Z^2 \cdot (-13.6 \text{ eV}) + \frac{5}{4} Z \cdot 13.6 \text{ eV}.
$$

## Variational Method Estimation

Lets try getting an even better estimation with the variational method. If we change the atom number $Z$ to a parameter $\alpha$ we get:

$$
    \psi_0 (r_1, r_2) = \frac{\alpha^{3}}{\sqrt{\pi} a_0^{3}} e^{\alpha (r_1+r_2)/a_0}.
$$

We now use this wave function as the base to calculate the expectation value of the Hamiltonian:

$$
    E(\alpha) = \langle \psi_0 (r_1, r_2) | \hat{H} | \psi_0 (r_1, r_2) \rangle = (-2 \alpha^2 + 4 Z \alpha - \frac{5}{4} \alpha) \cdot (-13.6 \text{ eV}).
$$

Here we can see that if we choose $\alpha = Z$ we get the same result as for the first order correction perturbation theory. The variational method says $E_\text{ground} \leq E(\alpha)$ so lets check where the minimum of $E(\alpha)$ is:

$$
    \frac{\partial E(\alpha)}{\partial \alpha} = (-4 \alpha + 4Z - \frac{5}{4}) = 0 \Longrightarrow \alpha = Z - \frac{5}{16}.
$$

So our estimated stationary state energy with the variational method is

$$
    E_\text{ground} = (Z-\frac{5}{16})^2 \cdot (-13.6 \text{ eV}).
$$

For example, the estimated stationary state energy with the variational method for helium ($Z=2$) is $-77.5$ eV which is close to the experimental value of $-78.975$ eV, in other words this method of estimation works very well! You now know a couple of methods to estimate the stationary state energy for electrons!

FILTERED POSTS

Allowed, Forbidden, Fermi, and Gamow-Teller Beta Decay Transition

Precession Rate of a Rolling Disk

Spherical Symmetric Charge Distribution of Nucleis

The ARMA Process

S-parameters in Circuits and Antennas

RECENT POSTS