Chakradhar Rangi

Example Talk

Sat, 01 Jun 2030 13:00:00 +0000

Click on the Slides button above to view the built-in slides feature.

Slides can be added in a few ways:

Create slides using Wowchemy’s Slides feature and link using slides parameter in the front matter of the talk file
Upload an existing slide deck to static/ and link using url_slides parameter in the front matter of the talk file
Embed your slides (e.g. Google Slides) or presentation video on this page using shortcodes.

Further event details, including page elements such as image galleries, can be added to the body of this page.

Classical Benchmarks of a Symmetry-Adapted Variational Quantum Eigensolver for Real-Time Green's Functions in Dynamical Mean-Field Theory

Tue, 03 Feb 2026 00:00:00 +0000

Real-Time Iteration Scheme for Dynamical Mean-Field Theory: A Framework for Near-Term Quantum Simulation

Tue, 27 Jan 2026 00:00:00 +0000

Beyond Black and White: Adapting Models to Visual Domain Shift

Thu, 04 Dec 2025 00:00:00 +0000

This project explores Unsupervised Domain Adaptation (UDA) techniques to address covariate shift in neural networks. We demonstrate how models trained on a source domain (grayscale MNIST digits) fail to generalize to a visually complex target domain (colored and textured MNIST-M digits), even when the underlying classification task is identical.

The project is implemented across three progressive stages:

Baseline: Demonstrating the impact of domain shift on model performance.
Distance-based Alignment: Using Maximum Mean Discrepancy (MMD) to align latent feature distributions.
Adversarial Training (DANN): Implementing Domain-Adversarial Neural Networks with a Gradient Reversal Layer (GRL) for superior feature alignment.

The Problem: Covariate Shift

Neural networks typically assume that training and deployment data follow the same distribution. In this project, we utilize:

Source Domain : MNIST (grayscale digit images).
Target Domain : MNIST-M (digits blended with random color patches from the BSDS500 dataset).

While the label space remains the same, the visual shift causes a standard CNN to struggle with learning dataset-invariant features.

Methodology & Results

Task 1: Baseline (Domain Shift Demonstration)

A standard CNN was trained exclusively on the source domain (MNIST). While it achieved high accuracy on the source test set, performance plummeted by over 50% when evaluated on the unseen MNIST-M target domain, proving the model’s inability to generalize under covariate shift.

Task 2: Distance-based Domain Adaptation

We implemented a “tug-of-war” optimization using Energy Distance (mathematically equivalent to MMD with a Gaussian kernel).

Approach: Forced the network to output latent representations that are statistically identical for both domains.
Results: Performance on MNIST-M jumped from ~43% to ~63%.
Key Finding: Statistical alignment is highly sensitive to batch size. Larger batch sizes (up to 1024) significantly improved statistical distribution estimates, peaking at 79.36% accuracy with tuned learning rates.

Task 3: Adversarial Domain Adaptation (DANN)

Inspired by GANs, we implemented a Domain-Adversarial Neural Network (DANN).

Architecture: Modified the CNN to include a Domain Discriminator and a Gradient Reversal Layer (GRL).
Mechanism: The feature extractor actively learns to “confuse” the discriminator, ensuring features are discriminative for classification but uninformative regarding the domain.
Results: This method achieved the highest robustness, with a mean accuracy of 89.04% using a fixed alpha strategy.

Results summary

Key Technical Concepts

Gradient Reversal Layer (GRL): Acts as an identity function during the forward pass but multiplies gradients by during backpropagation to facilitate adversarial training.
Unsupervised Adaptation: No target labels were ever used during the training process; the model adapts purely through feature alignment.

Quantifying Negative Convexity in Mortgage-Backed Securities

Sat, 15 Nov 2025 00:00:00 +0000

This project develops a quantitative framework to model the cash flows and price sensitivity of Mortgage-Backed Securities (MBS). The primary focus is investigating the phenomenon of Negative Convexity—the asymmetric risk profile caused by the embedded homeowner prepayment option.

Key Features

Custom Prepayment Model: Built a Python-based simulation engine for MBS valuation, featuring a path-dependent S-Curve model to quantify the impact of interest rate volatility and prepayment speeds on portfolio convexity.
Comparative Benchmarking: Quantified the “Convexity Gap” by comparing MBS performance against a synthetic non-callable bond.
Scenario Stress Testing: Analyzed portfolio Present Value (PV) sensitivity under ±25 and ±50 basis point interest rate shocks.

The Quantitative Challenge: Negative Convexity

Unlike standard fixed-income instruments, MBS investors are effectively Short Volatility. As rates drop, prepayments spike (prepayment risk), and as rates rise, prepayments dry up (extension risk).

Mathematical Foundation: The Geometry of Risk

To quantify how the Mortgage-Backed Security price ($P$) responds to interest rate fluctuations ($\Delta y$), we utilize a second-order Taylor expansion. This allows us to decompose price sensitivity into linear and non-linear components.

$$\frac{\Delta P}{P} \approx -D_{mod} \cdot \Delta y + \frac{1}{2} C \cdot (\Delta y)^2$$

1. First-Order Sensitivity: Modified Duration ($D_{mod}$)

Modified Duration represents the linear sensitivity of the bond price to a change in yields. In physics terms, this is the normalized first derivative—effectively the “velocity” of price change:

$$D_{mod} = -\frac{1}{P} \frac{\partial P}{\partial y}$$

2. Second-Order Sensitivity: Convexity ($C$)

Convexity captures the “curvature” of the price-yield relationship. It is the normalized second derivative—effectively the “acceleration” of price change:

$$C = \frac{1}{P} \frac{\partial^2 P}{\partial y^2}$$

3. The Negative Convexity Regime

For standard fixed-income instruments (e.g., non-callable Treasuries), $C > 0$. This Positive Convexity is a desirable trait, as it implies that prices rise faster than they fall for a given basis point move.

In the case of Mortgage-Backed Securities, the embedded homeowner call option introduces Negative Convexity ($C < 0$). As market yields decrease, the probability of prepayments spikes, causing principal to be returned prematurely. This results in a concave price-yield curve where:

$$\frac{\partial^2 P}{\partial y^2} < 0$$

This project benchmarks the MBS against a synthetic non-callable “Bullet” bond to visualize and quantify this Convexity Gap, measuring the cost of the “short volatility” position inherent in the mortgage market.

Results Summary

Our analysis demonstrates that the MBS exhibits significant “price capping” in falling-rate environments. At a -50bps shock, the MBS portfolio underperformed the benchmark by 3.83%, illustrating the substantial cost of the embedded call option.

Technical Skills Demonstrated

Language: Python (NumPy, Pandas, Matplotlib)
Finance: Fixed Income, Convexity/Gamma, Prepayment Modeling, Yield Curve Shifting.

Bee Health Early Warning System - Erdos Institute Data Science Bootcamp

Wed, 05 Nov 2025 00:00:00 +0000

Motivation and Objective

Honey bee colonies, which are essential for pollinating over $15 billion in U.S. crops annually, are collapsing at alarming rates. This poses a direct threat to national food security and economic stability. Currently, data on colony loss is a lagging indicator, meaning stakeholders can only react after a collapse has already occurred. In Phase 1 of our project, we addressed this by attempting to build a macro-level, “top-down” early warning system. The objective was to predict high or low colony loss at a U.S. state-by-state, quarterly level. However, this approach faced a significant challenge: averaging weather and pathogen data over vast geographic areas and entire quarters diluted the critical, localized patterns that truly stress a colony. Our models struggled to find a clear signal in this noisy data, which was reflected in poor F1 and ROC-AUC scores.

The challenges from Phase 1 directly motivated our pivot to Phase 2: a “bottom-up” approach to solve the problem at the scale where interventions actually happen - the individual hive. Our new objective is to create a model that predicts hive health (healthy/unhealthy) using highly localized weather data from the past one to two weeks. This approach is more robust for two key reasons: first, it aligns directly with the granular, local data that is available, and second, it allows us to use a standardized definition of “hive health” as established in recent research (Lower et al., 2024).

Stakeholders

This diagnostic model provides direct, tactical value to the people who can most immediately protect bee populations. Our primary stakeholders are Commercial Beekeepers, who can use this model to monitor their operations, enabling them to inspect and treat at-risk hives immediately rather than discovering a loss later. This tool would also be invaluable to Beekeeping Technology Companies (e.g., “Smart Hive” developers), as our model can serve as the software brain for their hardware sensors, turning raw data into an actionable warning system. Finally, the model would benefit Hobbyist Beekeepers by providing an accessible tool to help them manage hive health, as well as Pollination Service Providers who rely on a stable and healthy supply of bees to meet agricultural contracts.

Key Performance Indicators

To evaluate our classification models, we focused on three key performance indicators (KPIs): Accuracy, ROC-AUC, and the F1-score. We used these metrics to compare our more advanced models against each other and against a baseline model to quantify our improvements. While Accuracy provides a general measure of correctness, it can be misleading in a dataset where ‘healthy’ hives are less frequent but critical class to identify. Therefore, the ROC-AUC score was essential, as it measures a model’s fundamental ability to distinguish between the ‘healthy’ and ‘unhealthy’ classes. Most importantly, the F1-score was our primary metric for real-world performance, as it balances the critical trade-off between Precision (minimizing false alarms that waste a beekeeper’s time) and Recall (our main objective: successfully finding as many truly ‘unhealthy’ hives as possible).

Phase 1:

Data Sources and Processing

We created a comprehensive dataset for analyzing bee health by integrating three sources from 2015-2025. First, bee colony loss data (including specific causes like Varroa mites) from USDA NASS was aggregated by state and county. Second, this was merged with quarterly-averaged weather data (e.g., precipitation, temperature) from NOAA. Third, we added pathogen prevalence data, including Varroa counts and key virus levels, from USDA APHIS. This final, harmonized dataset aligns bee colony losses with corresponding weather and pathogen conditions by state, county, and quarter, providing a robust foundation for statistical analysis and predictive modeling.

Methodology / Modeling

Our Phase 1 objective was to predict future colony loss at the state level. To do this, we engineered a target variable by converting the raw quarterly loss percentage into a binary classification (High/Low Loss) based on the historical median loss for that quarter. This created a balanced target suitable for classification modeling. Our validation strategy was critical. Because the data was a time series, we used TimeSeriesSplit cross-validation. This ensured our models were always trained on past data (e.g., 2015-2020) to predict future outcomes (e.g., 2021), preventing data leakage. We experimented with several models, including LogisticRegression, XGBClassifier, and LGBMClassifier, to see if a signal could be found. We also conducted rigorous feature engineering, including lagging features and testing various preprocessing pipelines (imputation, scaling) to ensure a fair comparison.

Results and Conclusions

Our models conclusively demonstrated that predicting future state-level loss is not feasible with this data. Despite rigorous testing, our predictive models (forecasting Y(Q2) from X(Q1)) failed to find a signal, achieving an average ROC-AUC score of approximately 0.53—statistically indistinguishable from a random guess. Our critical analysis revealed four root causes for this failure:

Rapid Signal Decay: Our investigation showed that while a weak contemporaneous signal exists (linking Varroa in Q1 to loss in Q1), this signal decays too rapidly. The one-quarter lag was too long, and the health status from one quarter had no predictive power on the next quarter.
Extreme Feature Sparsity: The public pathogen data, which we hypothesized was a key predictor, was too sparse. Our analysis showed an average of only 4-8 samples per state, per quarter. This is not enough data to create a stable, representative feature, and our models correctly learned to treat it as statistical noise.
Destructive Aggregation: As stated in our objective, averaging weather and pathogen data over an entire state (like Texas) and a full quarter “dilutes” the signal. Localized events, which are the true drivers of colony health, are lost in this macro-level view.
Insufficient Data Depth: With only 40 quarters of data across ~50 states (roughly 2,000 total data points), the dataset was too shallow for a model to learn the weak, complex patterns hidden in the noise. This limited sample size made it impossible to train a generalizable model. This definitive failure was the primary motivation for our pivot to Phase 2. We concluded that a successful model must be built from the “bottom-up”, using granular, hive-level health data linked directly to local weather, which is precisely what Phase 2 accomplishes.

Phase 2:

Data Sources and Processing

Our model was built by combining two distinct datasets. First, we used Honeybee Colony Health Data from the “Predicting Honeybee Health” research article (Lower et al., 2024), which provided detailed hive inspection records (HCC_Inspections.csv) and apiary locations. Second, we used Daily Weather Summaries from the NOAA NCEI data portal for the corresponding apiary locations in North Carolina and Utah. To prepare the data for modeling, we converted the “Healthy” target variable to binary (1/0) and merged the two sources. Our primary feature engineering task was to prevent data leakage by creating time-series “lag” features (e.g., Prev_Health_Status) to capture the hive’s state from its previous inspection, as well as 7-day trailing weather features (e.g., Avg_tmax, Num_frost_days) to represent the environmental conditions leading up to the current inspection.

Methodology / Modeling

The nature of our data, which contains multiple inspections for the same hive over time, means that individual data points are not independent. A standard random split would cause data leakage, as the model could “memorize” a hive’s history from the training set. To simulate a real-world scenario, predicting health for a new, unseen hive, we used a GroupShuffleSplit on HiveID. This ensures all inspections for a single hive are confined to either the training or test set. We established a baseline with LogisticRegression and then experimented with ensemble models (RandomForest, LGBMClassifier, XGBClassifier). After hyperparameter tuning with RandomizedSearchCV, we selected the XGBClassifier as our final model. It yielded the strongest cross-validation ROC-AUC score (0.75), confirming it had the best ability to distinguish between ‘healthy’ and ‘unhealthy’ hives.

Results and Conclusions

Our final tuned XGBClassifier generalized well to the unseen holdout data, achieving a Test ROC-AUC of 0.78 and an Accuracy of 0.72. This strong predictive signal significantly outperformed our LogisticRegression baseline, which only achieved a Test ROC-AUC of 0.65. The model’s real-world value was confirmed by its Macro Average F1-score of 0.71, demonstrating a robust and balanced ability to identify both ‘healthy’ and ‘unhealthy’ hives. Feature importance analysis revealed that the model’s predictions were primarily driven by the hive’s recent history. The top three most important predictors were the Previous Queen Status (Prev_Queen_Status), whether it was the First Inspection (Is_First_Inspection), and the Previous Stressor Status (Prev_Stressors_Status). Key environmental factors, specifically the Number of Frost Days (Num_frost_days) and 7-day Average Temperature (Avg_tavg), also proved to be significant predictors of hive health.

Limitations and Future Work

While our Phase 2 model was successful (Test ROC-AUC: 0.78), its practical application is currently limited by several key factors. Our future work proposals directly address these limitations to build a more robust and scalable system.

Limitations

Limited Geographic and Climatic Generalizability: Our model was trained exclusively on data from North Carolina and Utah. These two regions, while different, do not represent the full spectrum of climatic and geographic conditions beekeepers face (e.g., the high humidity of Florida, the extreme cold of the upper Midwest, or the intensive migratory agriculture of California). The model would likely perform poorly in these new environments.
Limited Temporal Scope: The 2016-2019 data window is a narrow snapshot. The model has not been trained on major anomalous climate events, such as a widespread drought or an exceptionally harsh winter. Its ability to predict health outcomes during such extreme, non-linear events is untested and a significant risk.

Future Work

Expanded Feature Engineering: Our model used a 7-day trailing average for weather. This should be the subject of deeper experimentation. We hypothesize that different lag periods (e.g., 3-day for acute stress, 21-day for forage impact) could unlock new predictive signals.
Integrate Hive-Level Pathogen Data: Our Phase 1 project failed because state-level pathogen data was too sparse. The logical next step is to add hive-level pathogen screening (e.g., Varroa counts, DWV swabs) to our Phase 2 model. We predict these features would be immensely powerful, likely surpassing weather as a key secondary driver.
Transition from Manual to Automated Features: The ultimate goal is to move from a diagnostic tool to a real-time one. Future work should focus on replacing lagged features (Prev_Brood_Status) with data from in-hive sensors (e.g., temperature, humidity, acoustics, CO2). This would create a fully automated “early warning system” and is the natural path to a commercial product.
Scaling: This hive-level model acts as a precursor for solving the national-level problem. It proves that a “bottom-up" approach works. By deploying this model across more states, we can aggregate the predictions from thousands of individual hives. This would finally provide policymakers with a real, granular, and predictive tool to monitor national bee health, a goal our Phase 1 model could not achieve due to its “top-down” data limitations.

References

Lower, E., Kollaparthi, S. P., Rogers, R., Hassler, E., & Cazier, J. (2024, Jul 30). Predicting Honeybee Health: The Healthy Colony Checklist,Hive Scale and Weather Data. Data & Analytics for Good, (2). https://data-for-good.pubpub.org/pub/rg3364dl

Are financial stock returns normally distributed?

Sun, 24 Aug 2025 00:00:00 +0000

Have you ever heard that stock market returns (logarithmic returns to be precise) are supposed to follow that classic bell curve—the “normal distribution”? It’s one of the key assumptions of financial math, making a lot of complex theories work neatly. But the real world market data feels a lot messier than a perfect bell curve, doesn’t it?

So, in this blog, I decided to put this textbook theory to the test. Let’s grab some real stock data and see if it actually behaves the way the models say it should. We’ll be detectives, using a few simple graphical and statistical tools to uncover the truth.

Graphical tests

Visual checks are a great first step to get a feel for the data. Today, we will be looking at two such tests.

Quantile-Quantile(Q-Q) plot: With a Q-Q plot, the observed data’s quantiles are plotted against the theoretical quantiles (normal distribution in our case). This one sounds complicated, but it’s basically a “line-up” test. We ask our data’s values to line up against where they should be in a perfect normal distribution (theoretical quantiles). If our data is truly “normal,” the points on the plot will form a perfectly straight diagonal line. If they start wandering off—especially at the ends—it’s a dead giveaway that something else is going on.
Box plot: The box plot is a quick check for our data’s symmetry. Imagine a perfectly balanced scale: the median line is the center point, and the box and whiskers are evenly spread on both sides. If one side is much longer than the other or the median is off-center, it’s a clear signal that our data is skewed and not following a normal distribution.

The main drawback of these tests, however, is that they’re subjective—we’re sometimes just eyeballing it. They don’t give us a hard number to measure exactly how much our data deviates from normality. That’s where statistical tests come in - they give us a clear, quantifiable result.

Statistical tests

Statistical tests give us a clear mathematical verdict. They work by first assuming the data is normal and then calculating how likely it is we would see our particular dataset under that assumption.

This likelihood is called the p-value. In simple terms: if the p-value is low (less than 0.05), it means our data is a very bad fit for the normal distribution model. So, a low p-value gives us statistical proof that our data is not normal.

The specific test we’ll be using in this blog is the Shapiro-Wilk test. Another popular statistical test is the Kolmogorov-Smirnov test which we’ll not be using today.

Testing the evidence of normal distribution for JPMorgan Chase & Co

For this test, I’m choosing JPMorgan Chase & Co (JPM), but you can follow along with any stock you’re curious about! Just change the ticker from ‘JPM’ to something else (like ‘AAPL’ or ‘TSLA’) in the code below and see what you find.

Let us import the necessary python modules

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as stats
import seaborn as sns
import yfinance as yf 
import datetime as dt

# Setting global matplotlib parameters
import matplotlib as mpl
import matplotlib.dates as mdates 
mpl.rcParams["xtick.labelsize"] = 22
mpl.rcParams["ytick.labelsize"] = 22
mpl.rcParams['axes.titlesize'] = 22 
mpl.rcParams['figure.dpi'] = 300
mpl.rcParams['axes.labelsize'] = 24
mpl.rcParams['font.family'] = 'serif'
mpl.rcParams["mathtext.fontset"] = 'cm'

Downloading 10 year of JPM data

EndDate = dt.datetime.today()
StartDate = EndDate - pd.DateOffset(years=10)

JPMStockData = yf.download('JPM', start = StartDate, end = EndDate)["Close"]
JPMStockReturns = JPMStockData/JPMStockData.shift(1)
JPMStockLogReturns = np.log(JPMStockReturns.dropna())['JPM'].values
Dates = JPMStockReturns.dropna().index

Let us plot the Q-Q and box plot of log returns of JPM stock using built-in plotting functions from statsmodel and matplotlib

fig, axes = plt.subplots(1,2,figsize=(16, 6),dpi=300)

stats.probplot(JPMStockLogReturns, dist="norm", plot=axes[0])
axes[0].set_title("JPMorgan Chase & Co Log Returns Q-Q Plot")
axes[0].grid(True)
axes[1].boxplot(JPMStockLogReturns)
axes[1].set_xlabel("Variable")
axes[1].set_ylabel("Value")
axes[1].set_title("JPMorgan Chase & Co Log Returns Box Plot")
plt.tight_layout()
plt.show()

The Q-Q plot above shows the dots veering off the straight line, creating what analysts call “heavy tails.” Think of these as the market’s mood swings—the unexpectedly huge up and down days that the neat bell curve just can’t predict. Let us further perform statistical tests.

Shapiro-Wilk test

pvalue = stats.shapiro(JPMStockLogReturns)[1]
print("Normality Test Results of JPMStockLogReturns stock (Shapiro-Wilk test):")
print(f"p-value = {pvalue}")
if pvalue < 0.05:
    print("Statistically significant evidence that the data is NOT normally distributed.")
else:
    print("No statistically significant evidence against normality.")

Normality Test Results of JPMStockLogReturns stock (Shapiro-Wilk test):
p-value = 2.453415772115319e-39
Statistically significant evidence that the data is NOT normally distributed.

Our first verdict is here: Over the last 10 years, JPM’s returns are anything but normal, and both our statistical and graphical tests confirms it. The culprit? Our 10-year data captures wild market swings from events like COVID and global conflicts, leading to extreme data points. This sparks an interesting question: could there be a calmer period, free from extreme market volatility, where JPM log returns actually align with a normal distribution? To check for such periods, let us plot JPM log returns over the last 10 years.

fig,ax = plt.subplots(figsize=(10,6),dpi=300)
ax.plot(Dates,JPMStockLogReturns)
plt.title("Log returns of JPMorgan Chase & Co stock over a 10 year period")
ax.xaxis.set_major_locator(plt.MaxNLocator(6)) 
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
StartDate = dt.datetime(2021,1,1)
EndDate = dt.datetime(2022,1,1)
ax.axvspan(StartDate, EndDate, alpha=0.5, color='green')
plt.tight_layout()
plt.show()

Just look at this plot of the returns over time. You can literally see the wild swings from the COVID crash in 2020 and recent global conflicts.

But wait… see that green highlighted area in 2021? It looks much calmer. This got me thinking: maybe the normality assumption isn’t completely wrong, it just doesn’t work when the market is having a meltdown. What if we zoom in on that quiet period? Let’s investigate!

Let us download the data for JPMorgan Chase & Co during the highlighted area (2021)

JPMStockDataCustom = yf.download('JPM', start = StartDate, end = EndDate)["Close"]
JPMStockReturnsCustom = JPMStockDataCustom/JPMStockDataCustom.shift(1)
JPMStockLogReturnsCustom = np.log(JPMStockReturnsCustom.dropna())['JPM'].values
DatesCustom = JPMStockReturnsCustom.dropna().index

fig,ax = plt.subplots(figsize=(10,6),dpi=300)
ax.plot(DatesCustom,JPMStockLogReturnsCustom)
ax.xaxis.set_major_locator(plt.MaxNLocator(7)) 
ax.xaxis.set_major_formatter(mdates.DateFormatter('%y-%m'))
plt.title("Log returns of JPMorgan Chase & Co stock in 2021")
plt.tight_layout()
plt.show()

Q-Q and box plot - log returns of JPMorgan Chase & Co (2021)

fig, axes = plt.subplots(1,2,figsize=(16, 6),dpi=300)

stats.probplot(JPMStockLogReturnsCustom, dist="norm", plot=axes[0])
axes[0].set_title("JPM Log Returns Q-Q Plot (2021)")
axes[0].grid(True)
axes[1].boxplot(JPMStockLogReturnsCustom)
axes[1].set_xlabel("Variable")
axes[1].set_ylabel("Value")
axes[1].set_title("JPM Log Returns Box Plot (2021)")
plt.tight_layout()
plt.show()

Shapiro-Wilk test - Log returns of JPMorgan Chase & Co (2021)

pvalue = stats.shapiro(JPMStockLogReturnsCustom)[1]
print("Normality Test Results of JPMorgan Chase & Co stock (Shapiro-Wilk test):")
print(f"p-value = {pvalue}")
if pvalue < 0.05:
    print("Statistically significant evidence that the data is NOT normally distributed.")
else:
    print("No statistically significant evidence against normality.")

Normality Test Results of JPMorgan Chase & Co stock (Shapiro-Wilk test):
p-value = 0.7530741691589355
No statistically significant evidence against normality.

We are in for suprising results! The dots in the Q-Q plot have lined up and the p-value is 0.75. That’s a huge number, and it means we have no evidence to say the data isn’t normal. It seems we found a ’normal’ year for JPM. This is fascinating! It suggests that the normality assumption might hold up, but only in very specific, stable market conditions. This raises a compelling question for future exploration: what market conditions in 2021 fostered this normality, and can they be replicated in other stable periods? These questions are interesting future endeavors.

Further investigation: Trimming extremal values from the dataset

This got me thinking: we know the extreme market events are skewing our results. So what if we just removed them? If we trim off the most extreme daily returns, will the remaining “calmer” data finally follow a normal distribution?

To find out, we’ll use the box plot itself to define what counts as “extreme.” The plot’s “whiskers” create a boundary for the bulk of the data. Any data point that falls outside those whiskers is considered an outlier. We’ll chop off those pesky outliers and test the data that’s left.

We are re-plotting the JPM log returns over last 10 years to extract whiskers

fig, axes = plt.subplots(1,2,figsize=(16, 6),dpi=300)

stats.probplot(JPMStockLogReturns, dist="norm", plot=axes[0])
axes[0].set_title("JPMorgan Chase & Co Log Returns Q-Q Plot")
axes[0].grid(True)
WhiskersBoxPlot = axes[1].boxplot(JPMStockLogReturns)["whiskers"]
axes[1].set_xlabel("Variable")
axes[1].set_ylabel("Value")
axes[1].set_title("JPMorgan Chase & Co Log Returns Box Plot")
plt.tight_layout()
plt.show()

We are extracting the whiskers (threshold) from the box plot

LowerWhisker = WhiskersBoxPlot[0].get_ydata()[1]
UpperWhisker = WhiskersBoxPlot[1].get_ydata()[1]
print(f"The lower and upper threshold for extremal values are {LowerWhisker:.4f} and {UpperWhisker:.4f}")

The lower and upper threshold for extremal values are -0.0310 and 0.0326

Let us remove the extremal value from the log returns

Mask = np.where((JPMStockLogReturns <= UpperWhisker) & (JPMStockLogReturns >= LowerWhisker))
JPMStockLogReturnsTrimmed = JPMStockLogReturns[Mask]

Let us re-plot the JPM log returns over last 10 years omitting the extremal data

fig, axes = plt.subplots(1,2,figsize=(16, 6),dpi=300)

stats.probplot(JPMStockLogReturnsTrimmed, dist="norm", plot=axes[0])
axes[0].set_title("JPMorgan Chase & Co Log Returns Q-Q Plot")
axes[0].grid(True)
WhiskersBoxPlot = axes[1].boxplot(JPMStockLogReturnsTrimmed)["whiskers"]
axes[1].set_xlabel("Variable")
axes[1].set_ylabel("Value")
axes[1].set_title("JPMorgan Chase & Co Log Returns Box Plot")
plt.tight_layout()
plt.show()

Shapiro-Wilk test - Log returns of JPMorgan Chase & Co (Omitting the extremal values)

pvalue = stats.shapiro(JPMStockLogReturnsTrimmed)[1]
print("Normality Test Results of JPMorgan Chase & Co stock (Shapiro-Wilk test):")
print(f"p-value = {pvalue}")
if pvalue < 0.05:
    print("Statistically significant evidence that the data is NOT normally distributed.")
else:
    print("No statistically significant evidence against normality.")

Normality Test Results of JPMorgan Chase & Co stock (Shapiro-Wilk test):
p-value = 7.01930446211918e-07
Statistically significant evidence that the data is NOT normally distributed.

Interesting results! Even after omitting the extreme values, JPM’s log returns still refuse to conform to a normal distribution. This reveals a critical insight: the normality hypothesis does not just crumble under extreme market events—it falters even in calmer waters. We must critically analyze into factors like distribution skewness and question the assumption that daily log returns are independent and uncorrelated. While assuming normality simplifies the math, it strays far from reality, urging us to rethink our models to mimick the true stock behavior.

Conclusion

So, where does this leave us? My journey down this rabbit hole has shown me that the ’normal distribution’ of stock returns is more of a myth than a reality, at least over the long run. Even when we carefully removed the most extreme outliers, the data still refused to play by the rules. For me, the biggest takeaway is this: the outliers and the chaos in the data aren’t just noise to be ignored. They are the story. Relying on simplified models that assume a perfect, well-behaved world can be misleading.

Thanks for reading! I hope you enjoyed this blog as much as I did.

Interplay of non-Hermitian skin effect and electronic correlations in the non-Hermitian Hubbard model via Real-space dynamical mean field theory

Fri, 25 Jul 2025 00:00:00 +0000

Disorder Enhanced Thermalization in Interacting Many-Particle System

Sun, 25 May 2025 00:00:00 +0000

Real-Space Dynamical Mean-Field Theory (R-DMFT) for non-Hermitian strongly correlated systems

Sat, 08 Feb 2025 00:00:00 +0000

1. Background and Motivation

Non-Hermitian quantum systems have recently garnered significant interest due to their ability to describe open quantum systems, systems with gain and loss, or effective models for complex physical scenarios. Unlike traditional Hermitian systems, non-Hermitian Hamiltonians can exhibit complex eigenvalues and non-orthogonal eigenvectors, giving rise to fascinating phenomena such as the non-Hermitian skin effect (NHSE). The NHSE is characterized by the localization of a macroscopic number of eigenstates at the boundaries under open boundary conditions, challenging the conventional understanding of bulk-boundary correspondence.

The Hubbard model, a cornerstone of condensed matter physics, provides a powerful framework for studying strongly correlated electron systems, capturing phenomena such as Mott transitions, magnetism, and superconductivity. While traditionally Hermitian with symmetric hopping, the introduction of asymmetric hopping extends the Hubbard model into the non-Hermitian domain. This modification can represent biased particle flow, dissipation, or coupling to external environments, making it a versatile tool for exploring new physics.

We recently studied the non-Hermitian Hubbard model with asymmetric hopping, focusing on the interplay between the NHSE and strong electronic correlations. We employed the real-space dynamical mean-field theory (RDMFT) (also called inhomogeneous DMFT), a non-perturbative method renowned for its ability to accurately capture local correlation effects in strongly interacting systems. By analyzing the end-to-end Green’s function as probes of directional amplification, we demonstrate that strong correlations can suppress the skin effect, leading to a crossover from boundary-dominated to correlation-driven dynamics. Our systematic study reveals that correlations dominate at small to intermediate strength of asymmetric hopping, inducing exponential decay in the end-to-end Green’s function, but higher strength of asymmetric hopping can restore amplification even at strong interaction. These results illuminate the tunable interplay between correlations and non-Hermitian physics, suggesting avenues for engineering non-reciprocal transport in correlated open quantum systems.

However, here we highlight the R-DMFT algorithm for anyone interested!

2. Model and Methods

2.1 Non-Hermitian Hubbard Model

The Hamiltonian for the non-Hermitian Hubbard model on a 1D lattice is:

$$H = - \sum_{j,\sigma} \left( t_R c_{j+1,\sigma}^\dagger c_{j,\sigma} + t_L c_{j,\sigma}^\dagger c_{j+1,\sigma} \right) + U \sum_j n_{j\uparrow} n_{j\downarrow}$$

Where:

$c_{j,\sigma}^\dagger, c_{j,\sigma}$: Fermionic creation and annihilation operators at site $j$.
$\sigma = \uparrow, \downarrow$: Spin of the electron.
$t_R = t + \gamma$ and $t_L = t - \gamma$: Asymmetric hopping amplitudes, where $\gamma \neq 0$ denotes non-Hermiticity strength.
$U$: On-site Hubbard interaction strength.
$n_{j\sigma} = c_{j,\sigma}^\dagger c_{j,\sigma}$: Number operator.

The lattice has $L$ sites with open boundary conditions (OBC) to study the skin effect, where states localize at boundaries.

2.2 Real-Space DMFT Algorithm

This algorithm uses real-space DMFT to account for site-dependent properties, focusing on the retarded Green’s function $G^R(\omega)$ to compute the local density of states (LDOS). The self-energy is assumed local but varies across sites: $\Sigma_{ij}^R(\omega) = \delta_{ij} \Sigma_i^R(\omega)$.

Step 1: Initialize Parameters and Self-Energy
- Lattice size: Set $L$.
- Hopping parameters: Define $t_R = t + \gamma, t_L = t - \gamma$.
- Interaction: Set $U$.
- Frequency grid: Choose real frequencies $\omega \in [-\omega_{\text{max}}, \omega_{\text{max}}]$.
- Broadening factor: Set $\eta > 0$ (e.g., $\eta = 0.01$).
- Initial self-energy: For each site $i$, initialize $\Sigma_i^R(\omega)$ (e.g., $\Sigma_i^R(\omega) = 0$ or Hartree term $U n_{\text{opp}}$).
Step 2: Construct Non-Interacting Hamiltonian Build the non-interacting Hamiltonian $H_0$ as an $L \times L$ tridiagonal matrix for OBC: $H_0(i, i+1) = -t_R$ and $H_0(i+1, i) = -t_L$. Example for $L = 4$:

$$H_0 = \begin{pmatrix} 0 & -t_R & 0 & 0 \ -t_L & 0 & -t_R & 0 \ 0 & -t_L & 0 & -t_R \ 0 & 0 & -t_L & 0 \end{pmatrix}$$
Step 3: Compute the Lattice Green’s Function Form the diagonal self-energy matrix $\Sigma^R(\omega)$ and compute: $$G^R(\omega) = \left[ (\omega + i\eta) I - H_0 - \Sigma^R(\omega) \right]^{-1}$$
Step 4: Extract Local Green’s Functions For each site $i$, extract the diagonal element: $$G_{ii}^R(\omega) = \left[ G^R(\omega) \right]_{ii}$$
Step 5: Compute the Weiss Field and Hybridization Calculate the inverse Weiss field:

$$\mathcal{G}_{0,i}^R(\omega)^{-1} = G_{ii}^R(\omega)^{-1} + \Sigma_i^R(\omega)$$

The hybridization function is:

$$\Delta_i^R(\omega) = \omega - G_{ii}^R(\omega)^{-1} - \Sigma_i^R(\omega)$$
Step 6: Solve the Impurity Problem Solve the Anderson impurity model for each site. Using Iterated Perturbation Theory (IPT) at half-filling: $$\Sigma_i^R(\omega) = -i U^2 \int_{0}^\infty dt e^{i\omega t}[\beta_i(t)^2\alpha_i(-t) + \alpha_i(t)^2\beta_i(-t)]$$ where: $$\begin{Bmatrix} \alpha_i(t) \ \beta_i(t) \end{Bmatrix} = \int_{-\infty}^{\infty} \mathrm{d}\omega, e^{-i \omega t} \rho_{0,i}(\omega) f(\pm \omega)$$
Step 7: Update Self-Energy Mix the new and old self-energies for stability: $$\Sigma_i^R(\omega)^{\text{new}} = \alpha \Sigma_i^R(\omega)^{\text{computed}} + (1 - \alpha) \Sigma_i^R(\omega)^{\text{old}}$$
Step 8: Check for Convergence $$\text{error} = \max_{i,\omega} \left| \Sigma_i^R(\omega)^{\text{new}} - \Sigma_i^R(\omega)^{\text{old}} \right|$$ If $\text{error} < \text{tolerance}$, stop; otherwise, return to Step 3.
Step 9: Analyze the Skin Effect Compute the spectral function: $$A_i(\omega) = -\frac{1}{\pi} \Im G_{ii}^R(\omega)$$ The skin effect is observed if $A_i(\omega)$ shows enhanced values at boundary sites ($i=1$ or $i=L$).

3. Implementation of the Impurity Solver

As defined in the IPT approximation, the self-energy is computed through time-domain integrals. In numerical implementations, these Fourier transforms are carried out in discrete form using the Fast Fourier Transform (FFT), which scales as $\mathcal{O}(N\log N)$.

Discrete Fourier Transform (DFT):

$$\tilde{f}_k = \sum_{n=0}^{N-1}f_n e^{-i\frac{2\pi n k}{N}}$$

Inverse Discrete Fourier Transform (IDFT): $$f_n = \frac{1}{N}\sum_{k=0}^{N-1} \tilde{f}_ke^{i\frac{2\pi n k}{N}}$$

Out of Time Order Correlation of the Hubbard Model with Random Local Disorder

Tue, 23 Jul 2024 00:00:00 +0000

A measure for adiabatic contributions to quantum transitions

Fri, 22 Mar 2024 00:00:00 +0000

Engineering non-Hermitian Second Order Topological Insulator in Quasicrystals

Fri, 09 Feb 2024 00:00:00 +0000

crangi.github.io

Thu, 13 Jul 2023 00:00:00 +0000

A basic introduction to Majorana Edge modes in a Kitaev Chain

Sat, 06 Feb 2021 00:00:00 +0000

Welcome to another blog! This blog will give a very brief introduction and background of the Majorana zero modes and way to obtain them in a simple model known as the Kitaev Chain. The aim will be to translate the Kitaev Chain Hamiltonian into a Matrix form to obtain energy spectrum and edge modes for an open chain. We will obtain these Majorana zero modes at the edges of an open chain.

Majorana Fermions

In the year 1937, a new class of particles that are its own anti-particles were hypothesized by Ettore Majorana. This new class were termed Majorana fermions or particles for obvious reasons. They can be mathematically represented by creation and anhilation operators in second quantization:

$$\gamma^\dagger_k = \gamma_k$$ where, $\gamma_k$: anhilates a particle in quantum state k; $\gamma_k^\dagger$: creates a particle in quantum state k

The equation above conveys that the operators for creating and destroying the particles are same, which is the desired property of the Majorana in contrast to Dirac fermions for which they are distinct. The above operators shall be referred as Majorana operators. There are several conventions to define these Majorana operators. We will be adopting the following($f$ is the fermionic operator):

\begin{align} f^\dagger = (\gamma_1 + i\gamma_2)/2 \nonumber \qquad f = (\gamma_1 - i\gamma_2)/2 \label{2} \end{align}

Using this definition, the following identities can be checked

$\gamma_k^2 = 1$
$\left[\gamma_i,\gamma_j\right] = 2\delta_{ij}$

There have been lot of attempts in experimental high energy physics to hunt these particles.

Majorana Zero Modes

Let’s shift gears and enter the realm of condensed matter physics. In superconducting materials, certain quasiparticle(Bogoliubov) excitations display similar properties as that of Majorana fermions. The Majorana zero mode(MZM) is a special case of Bogoliubov quasiparticle operators which are linear combinations of electron and hole operator. The MZM satisfy the properties that they square to identity and are self-Hermitian.

If we look carefully at above equations, the two Majorana modes always occurs in pair for a single fermionic mode. Can we decouple them? This was cleverly answered by Professor Alexei Kitaev in 2001, through his Kitaev chain, which is a 1D p-wave superconducting model.

Kitaev Chain

Kitaev Chain is a tight binding lattice model of spinless electrons with nearest-neighbour hopping amplitude $t$, a p-wave superconducting pairing term $\Delta$ for nearest neighbours and on-site chemical potential $\mu$. The hamiltonian for a finite and open chain in second quantization can be jotted down \begin{equation} \hat{H} = \sum_{n=1}^{N-1}\left[t\left(f_n^\dagger f_{n+1} + f_{n+1}^\dagger f_n\right) + \Delta\left(f_n f_{n+1} + f_{n+1}^\dagger f_n^\dagger\right)\right] - \sum_{n=1}^{N} \mu \left( 2f_n^\dagger f_n - 1 \right) \label{real_space_Hamiltonian} \end{equation} where $t,\Delta$ and $\mu$ are all real and $\gamma > 0$ without loss of generality. Similar to the definition above, we define the following Majorana operators \begin{align} a_{2n-1} = f_n + f_n^\dagger \quad \text{and} \quad a_{2n} = i(f_n - f_n^\dagger) \label{4} \end{align} for n = 1,2,…,N. The Kitaev Chain Hamiltonian can be rewritten in terms of these Majorana operators \begin{equation} \hat{H} = i\sum_{n=1}^{N-1} \left[J_xa_{2n}a_{2n+1} - J_ya_{2n-1} a_{2n+2}\right] + i\sum_{n=1}^{N}\mu a_{2n-1}a_{2n} \end{equation} where, $J_x = \frac{1}{2}(t - \Delta) \quad \text{and} \quad J_y = \frac{1}{2}(t + \Delta)$

Let us carefully examine the Hamiltonian in the Majorana basis. A brief look will help us solve the problem of decoupling. Let us look at two special cases in the parameter space of $t, \Delta$ and $\mu$.

Case I: $t = \Delta$ and $\mu = 0$

The conditions clearly imply $J_x=0$ and $J_y=t$ which results in \begin{equation} \hat{H} = -i\sum_{n=1}^{N-1} \left[t a_{2n-1} a_{2n+2}\right] \label{6} \end{equation}

Did we end up finding something interesting? A close attention to the equation above reveals that the even(odd) Majorana mode at the left(right) end has vanished and it is decoupled from the rest of the modes. This is exactly what we set out for! The unpaired Majorana modes.

Case II: $t = -\Delta$ and $\mu = 0$

This time $J_y = 0$ and $J_x = t$ which results in \begin{equation} \hat{H} = i\sum_{n=1}^{N-1} \left[t a_{2n}a_{2n+1}\right] \end{equation}

This time, the odd(even) Majorana mode at the left(right) end of the chain is unpaired with the rest of the chain.

Wonderful! We know a way to obtain the unpaired Majorana end modes. Let us test this out through a Python Snippet!

Let us import the necessary Python modules

import numpy as np			# Numerical computation
import numpy.linalg as la		# Linear Algebra 
import matplotlib as mp			# Generating plots
import matplotlib.pyplot as plt
plt.style.use('seaborn')		# Setting the plotting style
mp.rcParams['figure.figsize'] = (9, 7)  # Setting the size of the plots
%matplotlib notebook			# Inline plotting

Let us define the various parameters of the system

Nsites = 25              # Number of lattice sites
Nprime = 2*Nsites
e_threshold = 1E-6	 # Threshold for finding zero eigenstates
params = {
't' : 2.0,               # Nearest neighbor hopping
'Delta' : 2.0,           # Superconducting pairing term
'mu' : 0.0               # Chemical potential
}

We define the Kitaev Chain Hamiltonian in the Majorana basis

def kitaev_ham(Nsites,params):
    Hmat = np.zeros([Nprime,Nprime])		# Declare a 2Nx2N matrix
    Jx = 0.5*(params['t'] - params['Delta'])	
    Jy = 0.5*(params['t'] + params['Delta'])
    
    for n in range(Nsites-1):
        Hmat[2*n,2*n+1]   = Jx
        Hmat[2*n+1,2*n]   = -Jx
        Hmat[2*n-1,2*n+2] = -Jy
        Hmat[2*n+2,2*n-1] = Jy
        Hmat[2*n-1,2*n]   = params['mu']
        Hmat[2*n,2*n-1]   = -params['mu']

    Hmat[2*(Nsites-1)-1,2*(Nsites-1)] = params['mu']
    Hmat[2*(Nsites-1),2*(Nsites-1)-1] = -params['mu']
    Hmat = 1j*Hmat

    return Hmat

Now that we have defined our Hamiltonian in a Matrix representation, let us visualise it.

Visualize the Hamiltonian Matrix

def visualise_mat(Hmat):
    plt.imshow(Hmat.imag) 	# The real part of the matrix is a zero matrix
    plt.colorbar()
    plt.show()

We have fixed the system parameters such that we realise case I. Let us see if we get any zero energy modes.

Diagonalise the Hamiltonian matrix and plot the spectrum

def plot_spectrum(Hmat):
    evals,evecs = la.eigh(Hmat)
    evals = evals.real
    plt.scatter(np.arange(len(evals)),evals)
    plt.title('Energy Spectrum of Chain with {} Sites'.format(Nsites))
    plt.show()
    
plot_spectrum(kitaev_ham(Nsites,params))

Hurray! We do see two energy modes at zero energy. Rest of the spectrum lies above and below the zero energy and also well separated. We can take it further and verify if the eigenvectors corresponding to these zero energy modes are localised at the edges of the chain.

Check for zero energy modes

# Extract the indices of energy modes close to zero
def check_zeromodes(evals):
    nzmodes = 0
    zmodes_ind = np.where(abs(evals) <= e_threshold)[0]
    return zmodes_ind,len(zmodes_ind)

Plot the probability distribution of the Majorana zero modes(edge modes)

def plot_zeromodes(evals,evecs,params):

    param_info = '\n'.join((
    r'$\mu=%.2f$' % (params['mu']),
    r'$t=%.2f$' % (params['t']),
    r'$\Delta=%.2f$' % (params['Delta'])))

    zmodes_ind,cnt_zmodes = check_zeromodes(evals)
    if cnt_zmodes > 0:
        fig,ax = plt.subplots(1,cnt_zmodes,figsize=(20, 10))
        fig.suptitle('Probability distribution of Zero modes',fontsize=20, fontweight='bold')
        for cnt in range(cnt_zmodes):
            ax1 = ax[cnt]
            ax1.plot(np.abs(evecs[:,zmodes_ind[cnt]])**2)
            ax1.set_title('Edge mode {}'.format(cnt+1),fontsize=20)
            ax1.set_xlabel('Site Number',fontsize=20)
            ax1.set_ylabel('$|\psi|^2$',fontsize=20)
            ax1.text(0.43, 0.95, param_info, transform=ax1.transAxes, fontsize=16,
        verticalalignment='top', bbox=dict(boxstyle="square",facecolor="white"))
            ax1.tick_params(axis='both', which='major', labelsize=16)
        #plt.savefig('Edge_modes_Kitaev.pdf')
        plt.show()
        
evals,evecs = la.eigh((kitaev_ham(Nsites,params)))
plot_zeromodes(evals,evecs,params)

Indeed, we see that wavefunction($|\psi|^2$ to be precise) corresponding to these zero energy modes are localised at the edges of the chain. We can also repeat the same exercise for Case II.

Wonderful! We know a way to obtain the unpaired Majorana end modes. The Majorana modes in the each case have zero energy and are localised at the end of the chain. Thus, it makes sense to coin the term, Majorana zero modes(MZMs). It is fascinating that the MZMs appear when tweaking the parameters! However, there might be some skepticism here as we obtained MZMs at a very special points in the parameter space. Let us vary the parameters of our system and check the fate of these zero energy modes.

var_mu = np.linspace(0,4,101)
var_energy = np.zeros([len(var_mu),Nprime])

for i in range(len(var_mu)):
    var_energy[i] = la.eigh(kitaev_ham(Nsites,params = {'t' : 2.0,'Delta' : 2.0, 'mu' : var_mu[i]}))[0]

plt.title("Energy Spectrum as a function of $\mu/t$ (N = 15)")
for i in range(Nprime):
    plt.plot(var_mu,var_energy[:,i])
plt.ylabel('Energy')
plt.xlabel('$\mu/t$')
plt.show()

It is a very interesting result. The zero energy modes exist even when we tweak the chemical potential. In other words, the MZMs don’t just appear by just fine tuning to a special parameter values but isolated MZMs stay protected as long as the bulk gap is finite. The reason behind this exciting observation is due to a concept known as Symmetry Protected Topological Phases(SPTP) which is a story for some other blog. For more details on this, refer to the following online course - Topology in Condensed Matter.

Thank you for your time! Happy learning!

Heat capacity of solids

Thu, 23 Jul 2020 00:00:00 +0000

In this parts of the notes, we will try to understand the historical evolution in trying to explain the heat capacity of solids. I am following ‘The Oxford Solid state Basics’ book by Prof. Steven H. Simon. If the material covered is unclear, please refer the book suggested above. In any case, go through that book. It is a lovely book to read!

Historical note: It was well known that the heat capacity per atom for a solid is given by $C \approx 3k_B$ or in other terms, molar heat capacity is $C \approx 3R$, where $R$ is the Universal gas constant. This empirical observation was known as Dulong-Petit Law.

Let us look at heat capacities for different chemicals at room temperature($T = 25^oC$)¹:

The above figure shows the heat capacities for various chemical elements at a temperature of 25 $^0C$ and a pressure of 100kPa. It is evident that most of them follow Dulong-Petit Law.

How do we explain this empirical observation? The answer was provided by Boltzmann using classical statistical mechanics. We will take his approach and understand the source of this empirical law.

Boltzmann approach: Classical Statistical Mechanics

The key idea in his approach was to assume atoms inside a solid to be in a harmonic potential under the influence of neighbouring atoms:

Imagine the same model with chain of such atoms forming a solid in 3 spatial dimensions(3D) at a fixed temperature T. Now, the Hamiltonian for each atom can be written as: $$ \mathcal{H}(\mathbf{q},\mathbf{p}) = \frac{\mathbf{p}^2}{2m} + \frac{1}{2}m\omega^2\mathbf{q}^2$$

Note:

We assumed the spring constant to be same in all three dimensions.

Quantities $\mathbf{p}$ and $\mathbf{q}$ are three dimensional vectors.

Since we have a large collection of such atoms at a fixed temperature T, we can find its energy and heat constant using classical canonical partition function. The single-particle partition function is given by, $$ Z_1 = \frac{1}{h^3}\int\int e^{-\beta \mathcal{H}(\mathbf{q},\mathbf{p})}d\mathbf{q}d\mathbf{p}$$

We shall explicitly evaluate this quantity. It will be handy to recall the Gaussian integral of one variable: $$\int_{-\infty}^{\infty} e^{-ax^2}dx = \sqrt{\frac{\pi}{a}}$$

Let us get our hands dirty!

$$Z_1 = \frac{1}{h^3}\int\int e^{-\beta(\frac{\mathbf{p}^2}{2m} + \frac{1}{2}m\omega^2\mathbf{q}^2)}d\mathbf{q}d\mathbf{p} \ = \frac{1}{h^3} \int exp{\frac{-\beta\mathbf{p}^2}{2m}} d\mathbf{p} \int exp{\frac{-\beta m\omega^2}{2}\mathbf{q}^2} d\mathbf{q} $$

We can further separate the 3D integral into a product of three 1D integrals in both momentum and position coordinates.

$$ \begin{align} Z_1 = \frac{1}{h^3} \int_{-\infty}^{\infty} exp{\frac{-\beta{p_x}^2}{2m}} dp_x \int_{-\infty}^{\infty} exp{\frac{-\beta{p_y}^2}{2m}} dp_y \int_{-\infty}^{\infty} exp{\frac{-\beta{p_z}^2}{2m}} dp_z \\\ \text{x} \int_{-\infty}^{\infty} exp{\frac{-\beta m\omega^2q_x^2}{2}} dq_x \int_{-\infty}^{\infty} exp{\frac{-\beta m\omega^2q_y^2}{2}} dq_y \int_{-\infty}^{\infty} exp{\frac{-\beta m\omega^2q_z^2}{2}} dq_z \end{align} $$

$$= \frac{1}{h^3} \left[\int_{-\infty}^{\infty} exp{\frac{-\beta{p_x}^2}{2m}} dp_x\right]^3 \left[\int_{-\infty}^{\infty} exp{\frac{-\beta m\omega^2q_x^2}{2}} dq_x\right]^3$$ Using the 1D Gaussian integrals above, we get

$$ Z_1 = \frac{1}{h^3} \left[\sqrt{\frac{2\pi m}{\beta}}\right]^3 \left[\sqrt{\frac{2\pi}{\beta m \omega^2}} \right]^3 = \left[\frac{k_B T}{\hbar\omega}\right]^3$$ where $k_B$ is the Boltzmann constant and $\hbar = h/2\pi$. The N particle partition function can be obtained by assuming that the particles are distinguishable: $$Z_N = [Z_1]^N = \left[\frac{k_B T}{\hbar\omega}\right]^{3N}$$

We know, from the connection to thermodynamic quantities, the internal energy of this system can be obtained by $$ E = -\frac{\partial}{\partial\beta}lnZ_N = 3Nk_BT$$

Thus, the heat capacity per atom is $$C_v = \frac{\partial}{\partial T} \frac{E}{N} = 3k_B$$

Hooray! We have obtained our result of $3k_B$ for the specific heat. Boltzmann’s idea worked and we can now explain this empirical observation. Thus, Boltzmann accounted lattice vibrations as the source for the heat capacity and modelled it classically as masses connected by a spring.

Is it the end of the story? Even the graph above will tell you that there are outliers to this empirical law. Particularly, heat capacity of diamond deviates from Dulong-Petit law at room temperatures. Also, heat capacity strongly depends on temperature, and becomes smaller at lower temperatures(Apparently room temperature is a “low” temperature for diamond). What are we missing here? Einstein came up with a brilliant idea to consider quantum mechanics as the key here.

Einstein’s approach: Quantum Harmonic Oscillator

Einstein followed a similar approach to that of Boltzmann except, he used quantized harmonic oscillator levels for each atom, given by $$\epsilon_n = \hbar \omega(n+3/2)$$

where $n = n_x+n_y+n_z$. The partition function in this case in given by: $$ \begin{align} Z_1 &= \sum_{n_x,n_y,n_z\geq0}e^{-\beta \epsilon_n} = \sum_{n_x,n_y,n_z\geq0}e^{-\beta \hbar \omega(n_x+n_y+n_z+3/2)} \\\ &= e^{-3\beta \hbar\omega/2}\sum_{n_x\geq0}e^{-\beta\hbar\omega n_x} \sum_{n_y\geq0}e^{-\beta\hbar\omega n_y} \sum_{n_z\geq0}e^{-\beta\hbar\omega n_z} \end{align} $$

Each sum represents a geometric series and thus reduces to $$ Z_1 = \left(\frac{e^{-\beta \hbar\omega/2}}{1-e^{\beta \hbar\omega}}\right) \left(\frac{e^{-\beta \hbar\omega/2}}{1-e^{\beta \hbar\omega}}\right)\left(\frac{e^{-\beta \hbar\omega/2}}{1-e^{\beta \hbar\omega}}\right) = \left[\frac{1}{2sinh(\beta\hbar\omega/2)}\right]^3$$

The N particle partitions function is thus given as $$Z_N = \left[\frac{1}{2sinh(\beta\hbar\omega/2)}\right]^{3N}$$

and energy $$ E = -\frac{\partial}{\partial\beta}lnZ_N = 3N\frac{\hbar\omega}{2}coth\left(\beta\hbar\omega/2\right) = 3N\hbar\omega\left(n_B(\beta\hbar\omega) + \frac{1}{2}\right) = 3NE_1$$

where $n_B(x) = \frac{1}{e^{x}-1}$ is the Bose occupation factor. The quantity $E_1$ stands for energy of a single harmonic oscillator in 1D and it can be interpreted as $n_B$ bosonic modes of energy $\hbar \omega$ are occupied.

The heat capacity per atom can be similarly obtained as above: $$C_v = \frac{\partial}{\partial T} \frac{E}{N} = 3k_B(\beta\hbar\omega)^2\frac{e^{\beta\hbar\omega}}{(e^{\beta\hbar\omega}-1)^2}$$

Obervation: The plot on the left tells us that lower energy modes have higher occupation than the higher energy modes. The plot on the right gives us the energy as a function of time. It shows us that as $K_BT\ll\hbar\omega$, the heat capacity, which is the slope, tends to zero.

This can be further verified by explicitly working out the heat capacity per atom: $$C_v = \frac{\partial}{\partial T} \frac{E}{N} = 3k_B(\beta\hbar\omega)^2\frac{e^{\beta\hbar\omega}}{(e^{\beta\hbar\omega}-1)^2}$$

The expression above might not be intuitive at first sight. We can check whether we obtain Dulong-Petit law in the high temperature limit: $K_BT \gg \hbar \omega$ $$ \frac{\hbar\omega}{K_BT} \ll 1 \implies \text{exp}(\frac{\hbar\omega}{K_BT}) \approx 1 \ \implies (e^{\beta\hbar\omega}-1) \approx \beta\hbar\omega \quad(\text{using Taylor expansion of exponential}) $$

Using these approximations, the high-temperature asymptotic behaviour of heat capacity reduces to $C = 3K_B$. We have managed to recover Dulong-Petit Law!

Similarly, we can obtain the low temperature limit: $K_BT \ll \hbar \omega$ $$ \frac{\hbar\omega}{K_BT} \gg 1 \implies \text{exp}(\frac{\hbar\omega}{K_BT}) \gg 1 \ \implies (e^{\beta\hbar\omega}-1) \approx \text{exp}(\frac{\hbar\omega}{K_BT}) $$ Using these approximations, the low-temperature asymptotic behaviour of heat capacity reduces to $$ C_v = 3k_B(\beta\hbar\omega)^2e^{-\beta\hbar\omega} $$

We see that heat capacity drops exponentially as we lower the temperature. All our findings, can be captured by plotting heat capacity as a function of temperature:

This exponential drop can be understood by looking at the energy spectrum of each atom. As we lower the temperature below the characteristic energy, i.e., $k_BT < \hbar\omega$, all atoms start occupying the ground state or low lying states according Bose distribution. Since the spectrum is discrete, the thermal energy is not enough to keep them in higher states and energy drop is higher(exponential).

Thus, the Einstein’s model qualitatively captures the low temperature behaviour of heat capacity in solids. Let us look at how it performs quantitatively and compares with a experimental data.

https://solidstate.quantumtinkerer.tudelft.nl/1_einstein_model/ ↩︎

Band Structure and Edge modes of a Graphene Nano Ribbon

Mon, 13 Jul 2020 00:00:00 +0000

Quick Summary

We are interested in the edge modes of the Graphene Nano Ribbon. A nano ribbon is a 2D structure in which the molecules extend infinitely in one of the directions and finite in the other. One more visualizations of the above is forming a chain of atoms(Periodic boundary conditions) in one of the directions and finite chain(Open boundary conditions) in the other direction.

Importing the required libraries

import numpy as np
import matplotlib as mp
from matplotlib import pyplot as plt
import ipywidgets as ipy
plt.style.use('seaborn')
mp.rcParams['figure.figsize'] = (9, 7)
%matplotlib notebook

Modules for building the model and computing the eigenenergies and values

def gzig_zag(N = 25, t = 1.0, kval = 0.5):
    Hmat = np.zeros([2*N,2*N], dtype = complex)
    for i in range(len(Hmat)-1):
        Hmat[i][i+1] = (1-i%2)*t*(1+np.exp(1j*kval)) + (i%2)*t #0.5*np.random.random()
        Hmat[i+1][i] = np.conj(Hmat[i][i+1])
        #Hmat[i][i] = np.random.random()
    return Hmat

def compute(Hmat):
    eigenenergies, eigenstate = np.linalg.eigh(Hmat)
    return eigenenergies, eigenstate

Energy spectrum as a function of ‘k’

N = 15; t = 1.0
k = np.linspace(-4,4,1001)
eig = np.zeros([len(k),2*N])
for _ in range(len(k)):
    eig[_] = compute(gzig_zag(N,t,k[_]))[0]
plt.figure()
plt.plot(k,eig)
plt.xlabel("$k$"); plt.ylabel("$E(k)$")
plt.title("Energy Spectrum for the Zig-Zag Graphene Structure")
plt.show()

Visualization of the Edge modes

vecs = np.zeros([2*N,2*N], dtype = complex)
vecs = compute(gzig_zag(N,t,-3))[1]
plt.figure()
plt.plot(range(2*N),vecs[:,N], label="Zero energy edge mode")
plt.plot(range(2*N),vecs[:,0], label="Ground state")
plt.title("Eigenvectors of the Zig-Zag Structure of Graphene")
plt.legend()
plt.show()

Checking the Robustness of the edge modes

def gzig_zag_r(N = 25, t = 1.0, kval = 0.5):
    Hmat = np.zeros([2*N,2*N], dtype = complex)
    for i in range(len(Hmat)-1):
        Hmat[i][i+1] = (1-i%2)*(t-np.random.random())*(1+np.exp(1j*kval)) + (i%2)*(t+np.random.random()) #0.5*np.random.random()
        Hmat[i+1][i] = np.conj(Hmat[i][i+1])
        #Hmat[i][i] = np.random.random()
    return Hmat

plt.figure()
N = 25; t = 1.0
k = np.linspace(-4,4,1001)
eig = np.zeros([len(k),2*N])
for _ in range(len(k)):
    eig[_] = compute(gzig_zag_r(N,t,k[_]))[0]
plt.plot(k,eig)
plt.xlabel("$k$"); plt.ylabel("$E(k)$")
plt.title("Energy Spectrum for the Zig-Zag Graphene Structure")
plt.show()

plt.figure()
vecs = np.zeros([2*N,2*N], dtype = complex)
vecs = compute(gzig_zag_r(N,t,-3))[1]
plt.plot(range(2*N),vecs[:,N-1], label="Zero energy edge mode")
plt.plot(range(2*N),vecs[:,0], label="Ground state")
plt.title("Eigenvectors of the Zig-Zag Structure of Graphene")
plt.legend()
plt.show()

Slides

Tue, 05 Feb 2019 00:00:00 +0000

Create slides in Markdown with Wowchemy

Wowchemy | Documentation

Features

Efficiently write slides in Markdown
3-in-1: Create, Present, and Publish your slides
Supports speaker notes
Mobile friendly slides

Controls

Next: Right Arrow or Space
Previous: Left Arrow
Start: Home
Finish: End
Overview: Esc
Speaker notes: S
Fullscreen: F
Zoom: Alt + Click
PDF Export

Code Highlighting

Inline code: variable

Code block:

porridge = "blueberry"
if porridge == "blueberry":
    print("Eating...")

Math

In-line math: $x + y = z$

Block math:

$$ f\left( x \right) = ;\frac{{2\left( {x + 4} \right)\left( {x - 4} \right)}}{{\left( {x + 4} \right)\left( {x + 1} \right)}} $$

Fragments

Make content appear incrementally

{{% fragment %}} One {{% /fragment %}}
{{% fragment %}} **Two** {{% /fragment %}}
{{% fragment %}} Three {{% /fragment %}}

Press Space to play!

One Two Three

A fragment can accept two optional parameters:

class: use a custom style (requires definition in custom CSS)
weight: sets the order in which a fragment appears

Speaker Notes

Add speaker notes to your presentation

{{% speaker_note %}}

- Only the speaker can read these notes
- Press `S` key to view
  {{% /speaker_note %}}

Press the S key to view the speaker notes!

Themes

black: Black background, white text, blue links (default)
white: White background, black text, blue links
league: Gray background, white text, blue links
beige: Beige background, dark text, brown links
sky: Blue background, thin dark text, blue links

night: Black background, thick white text, orange links
serif: Cappuccino background, gray text, brown links
simple: White background, black text, blue links
solarized: Cream-colored background, dark green text, blue links

Custom Slide

Customize the slide style and background

{{< slide background-image="/media/boards.jpg" >}}
{{< slide background-color="#0000FF" >}}
{{< slide class="my-style" >}}

Custom CSS Example

Let’s make headers navy colored.

Create assets/css/reveal_custom.css with:

.reveal section h1,
.reveal section h2,
.reveal section h3 {
  color: navy;
}

Questions?

Ask

Documentation

Mon, 01 Jan 0001 00:00:00 +0000

Erdos Institute Quant Finance Bootcamp

Mon, 01 Jan 0001 00:00:00 +0000

MS Thesis

Mon, 01 Jan 0001 00:00:00 +0000