# Magnetism in Superconductors

In this project, we will analyse data collected by Ana-Elena Tutueanu where we will look for a magnetic signal from neutrons sent onto the super-conductor $\text{La}_{2-x}\text{Sr}_x\text{CuO}_4$ where $x=0.07$ at different temperatures.

We expect to see two gaussian signals at $q=1 \pm \approx0.1$, and want to investigate at which temperature we lose the signal is it around $T_c=16K$.

## Basics

First we want to import our libraries.

import numpy as np
import matplotlib.pyplot as plt

from scipy import stats
from scipy.optimize import curve_fit, fsolve

from IPython.display import display, Latex


Before we load the data, we want to expand the file ranges to have a list of all the files in a dataset for each dataset.

temps = [1.5,4,7,10,20,40]
file_ranges = [
[18536,18538],
[18642,18644],
[18645,18647],
[18605,18607],
[18608,18610],
[[18543,18545],[18578,18580]],
]

expand = lambda limits: list(range(limits,limits+1))

# expand the file ranges.
expanded_files = []
for limits in file_ranges:
if isinstance(limits, int):
expanded_files.append(expand(limits))
else:
sub_expanded = []
for sublimits in limits:
sub_expanded.extend(expand(sublimits))

expanded_files.append(sub_expanded)


Next we define some settings and helper functions to load our files.

data_dir = "data"

filename = lambda number: f"{data_dir}/0{number}"

qcol = 2
ncol = 8

q = 1
delta = 0.1


To check that we load the files correctly, we plot all the files for each temperature.

fig,axes = plt.subplots(2,3, figsize=(14,8))

axes_flat = np.ndarray.flatten(axes)

for i,ax in enumerate(axes_flat):
ax.set_title(f"T={temps[i]}K")

ds = [loadfile(filename(expanded_files[i][j])) for j in range(len(expanded_files[i]))]

for j,data in enumerate(ds):
ax.plot(data[qcol],data[ncol],label=f"{expanded_files[i][j]}")

ax.set_xlabel("$q [q_0]$")
ax.set_ylabel("N over $122$s")
ax.legend()

axes_flat[-1].set_ylabel("N over $61$s")

fig.suptitle("Raw data") We see the $40K$ dataset has twice the number of files, but run for half the amount of time, so we need to create a combined $40K$ dataset.

## Investigate overlap in $T=40K$

It appears that $i$th file in each of the two ranges measure the same values. Which when combined gives you $3$ combined ranges. However, we want to know if we can assume the $i$th $q$ value are the same.

mean_error_q_interval = []
for i in range(3):
mean_error_q_interval.append(np.mean(np.abs(data1[qcol]- data2[qcol])))

for i,err in enumerate(mean_error_q_interval):
display(Latex("Average absolute difference between $q$ values in the {}'th subset {:.0E}.".format(i+1,err)))


Average absolute difference between $q$ values in the 1'th subset 9E-05

Average absolute difference between $q$ values in the 2'th subset 2E-05

Average absolute difference between $q$ values in the 3'th subset 1E-04

Since the absolute differences are so low, we will assume the $q$s match up, and we can now create a combined $40$K dataset by pair-wise adding the counts in the files of the two series.

for i in range(3):
file_numbers = [expanded_files[j] for j in [i,i+3]]

plt.plot(data1[qcol],data1[ncol]+data2[ncol], label=f"{file_numbers}+{file_numbers}")

plt.title(f"T=40K")
plt.xlabel("$q [q_0]$")
plt.ylabel("N over $122$s")
plt.legend()
plt.show() ## Final data

We can now create $6$ datasets - one for each of the measure temperatures - where we combine the three ranges of $q$ we have measured.

datasets = []

# handle first 5 datasets
for files in expanded_files[:-1]:
ds = []

for number in files:

data_combined = np.concatenate(ds,axis=1)

datasets.append(data_combined)

# handle 40K dataset
ds = []
for i in range(3):

data1[ncol] += data2[ncol]
ds.append(data1)

data_combined = np.concatenate(ds,axis=1)
datasets.append(data_combined)


### Plot raw data with combined datasets

fig,axes = plt.subplots(2,3, figsize=(14,8))
axes_flat = np.ndarray.flatten(axes)

for i,ax in enumerate(axes_flat):
ax.set_title(f"T={temps[i]}K")
ax.errorbar(datasets[i][qcol], datasets[i][ncol],yerr=np.sqrt(datasets[i][ncol]),fmt=".",label="Raw data")
ax.set_xlabel("$q [q_a]$")
ax.set_ylabel("N over $122$s")
ax.set_ylim(11000,16000)
ax.legend()

fig.suptitle("Combined datasets.") ## Subtract background and fit

We can now subtract the $40K$ background dataset from the other $5$ datasets. We note that the uncertainty of the $i$th datapoint for each temperature is $\sigma_i = \sqrt{N_i + N_{i,40K}}$ where $N_{i,40K}$ is the $i$th count of the $40K$ background dataset.

We fit the signals with two gaussians on top of a linear background. The gaussians are constrained such that they use the same $\sigma$ since this is dominated by the uncertainty of the equipment, and we constrain it further by forcing the centers $\mu=1\pm\delta$.

# fit function pieces
straight_line = lambda x, a,b: a*x+b
gauss = lambda x, A, mu, sigma: A * np.exp(-0.5*((x - mu) / sigma)**2.0)

# setup plots and axies
fig,axes = plt.subplots(2,3, figsize=(14,8))
axes = [plt.subplot2grid((2,6), pos, colspan=2)
for pos in [(0,0),(0,2),(0,4),(1,1),(1,3)]]

#      a    b    delta  sig   A    A
p0 = [[300, 200, 0.1,   0.04, 500, 500],
[300,-100, 0.1,   0.04,   0,   0],
[300,-300, 0.05,  0.04, 500, 500],
[300,-300, 0.05,  0.04, 200, 200],
[300,-300, 0.04,  0.02, 200, 200]]

popts = []
psigs = []

for i,ax in enumerate(axes):
ax.set_title(f"T={temps[i]}K")

qs = datasets[i][qcol]

signal_ns = datasets[i][ncol]
bg_ns = datasets[-1][ncol]
ns = signal_ns-bg_ns

# The new sigma is based on the combined
sigma_ns = np.sqrt(signal_ns+bg_ns)

ax.errorbar(qs, ns, yerr=sigma_ns, fmt=".",label="Data - Background")
ax.set_xlabel("$q [q_a]$")
ax.set_ylabel("N over 122s")
ax.set_ylim(-750,1600)

# linear fit
linpopt,linpcov = curve_fit(straight_line, qs, ns, p0=p0[i][:2], sigma=sigma_ns,maxfev=10000, absolute_sigma=True)
linpsig = np.diag(np.sqrt(linpcov))

# double gauss fit
model = lambda x, *p: straight_line(x,p,p) + \
gauss(x,p,1-p,p)  + \
gauss(x,p,1+p,p)

popt,pcov = curve_fit(model, qs, ns, p0=p0[i], sigma=sigma_ns,maxfev=10000, absolute_sigma=True)
psig = np.diag(np.sqrt(pcov))

popts.append(popt)
psigs.append(psig)

chi2r = lambda x,y,s,f,p: np.sum((y-f(x,*p))**2/s**2)/(len(x)-len(p))

chi2r_model = chi2r(qs,ns,sigma_ns,model,popt)
chi2r_line = chi2r(qs,ns,sigma_ns,straight_line,linpopt)

ax.text(0.73,1200,"$A_1$={:.0f}$\pm${:.0E}".format(popt[-2],psig[-2]))
ax.text(0.73,1050,"$A_2$={:.0f}$\pm${:.0E}".format(popt[-1],psig[-1]))
ax.text(0.73,800,"$\chi^2_m/n$={:.4f}".format(chi2r_model))
ax.text(0.73,650,"$\chi^2_l/n$={:.4f}".format(chi2r_line))

xs=np.linspace(min(qs), max(qs))

ax.plot(xs,straight_line(xs,*linpopt),'g', label="Linear fit",alpha=0.8)
ax.plot(xs,model(xs,*p0[i]),"--",c='orange',label="Initial Fit",alpha=0.45)
ax.plot(xs,model(xs,*popt),'orange', label="Model Fit")

ax.legend()


/opt/conda/envs/python3/lib/python3.7/site-packages/ipykernel_launcher.py:43: RuntimeWarning: invalid value encountered in sqrt
/opt/conda/envs/python3/lib/python3.7/site-packages/ipykernel_launcher.py:51: RuntimeWarning: invalid value encountered in sqrt We see that for the $1.5$K dataset, the two-gaussian model fits the data nicely. With the peaks being more than $2\sigma$ away from $0$, and the model $\chi^2_m/n$ is much close to $1$ than the linear fit $\chi^2_l/n$. We thus conclude that we do see a magnetic signal at $T=1.5K$.

If we look at the higher temperatures, it becomes more nuanced. We have greater $2\sigma$ peaks for both gaussians up to $7K$, and just around $2\sigma$ for $7K$ while just above $2\sigma$ again for $10K$. For $20K$, we clearly have less than $2\sigma$ peaks. The $\chi^2/n$ for the model and the line fits also become comparable somewhere between $10K$ and $20K$.

This would suggest that the critical temperature lies between $10K$ and $20K$.

## Investigations into the critical temperature $T_c$

Ideally, we would like to constrain our critical temperature to be more specific than somewhere between $10$ and $20K$. One way of doing this is by fitting the significance of the peaks $(A/\sigma)$ over temperature.

We the uncertainties for $A$ through curve_fit, and we can use Barlow 5.24 to estimate the error in $\sigma$

$\sigma_\sigma = \frac{\sigma}{\sqrt{2N}}$

Through error propagation, we find for $f(A,\sigma)=\frac{A}{\sigma}$

\begin{aligned}\sigma_f^2 &= \Big(\frac{df}{dA}\Big)^2 \sigma_A^2 + \Big(\frac{df}{d\sigma}\Big)^2 \sigma_\sigma^2 \\ &= \Big(\frac{1}{\sigma}\Big)^2 \sigma_A^2 + \Big(\frac{-A}{\sigma^2}\Big)^2 \sigma_\sigma^2 \end{aligned}

where $\sigma=\sigma_A$.

popts = np.array(popts)
psigs = np.array(psigs)

A1col, A2col = 4,5

fig, (ax1,ax2,ax3) = plt.subplots(1,3,figsize=(18,6))

ax1.plot(temps[:-1], popts[:,A1col]/psigs[:,A1col],".",c="tab:blue", label="$A_1/ \sigma_{A_1}$")
ax1.plot(temps[:-1], popts[:,A2col]/psigs[:,A2col],".",c="tab:orange", label="$A_2/ \sigma_{A_2}$")
ax1.set_title("Significance of $A$ without errorbars")
ax1.legend()

ax2.errorbar(temps[:-1], popts[:,A1col]/psigs[:,A1col],yerr=psigs[:,A1col], fmt=".", c="tab:blue", label="$A_1/ \sigma_{A_1}$")
ax3.errorbar(temps[:-1], popts[:,A2col]/psigs[:,A1col],yerr=psigs[:,A2col], fmt=".", c="tab:orange", label="$A_2/ \sigma_{A_2}$")
ax2.set_title("Significance of $A_1$ with errorbars.")
ax2.legend()
ax3.set_title("Significance of $A_2$ with errorbars.")
ax3.legend()

decay = lambda x, A, tau: A/(x**tau) #A*np.exp(-tau*x)

sigmas = lambda A, As, ss: np.sqrt((1/As)**2 *As**2 + (-A/As**2)**2 * ss**2)

sigmasA1,sigmasA2 = sigmas(popts[:,[A1col,A2col]],
psigs[:,[A1col,A2col]],
psigs[:,[A1col,A2col]]/np.sqrt(2*len(psigs[:,[A1col]]))).T

poptA1, pcovA1 = curve_fit(decay, temps[:-1], popts[:,A1col]/psigs[:,A1col], absolute_sigma=True, sigma=sigmasA1)
poptA2, pcovA2 = curve_fit(decay, temps[:-1], popts[:,A2col]/psigs[:,A2col], absolute_sigma=True, sigma=sigmasA2)

xs = np.linspace(1.5,20)
ax1.plot(xs,decay(xs,*poptA1),"tab:blue", label="fit $A_1/ \sigma_{A_1}$")
ax1.plot(xs,decay(xs,*poptA2),"tab:orange", label="fit $A_2/ \sigma_{A_2}$")

ax1.plot([1.5,20],[2,2], "r--", linewidth=1, alpha=0.5, label="$2\sigma$ confidence level")
ax2.plot([1.5,20],[2,2], "r--", linewidth=1, alpha=0.5, label="$2\sigma$ confidence level")
ax3.plot([1.5,20],[2,2], "r--", linewidth=1, alpha=0.5, label="$2\sigma$ confidence level")

ax1.set_xlabel("T")
ax1.set_ylabel("$A/ \sigma_{A}$")
ax2.set_xlabel("T")
ax2.set_ylabel("$A_1/ \sigma_{A_1}$")
ax2.set_xlabel("T")
ax2.set_ylabel("$A_2/ \sigma_{A_2}$")
plt.show() # when does A1 intersect the 2sigma confidence level?
x0 = fsolve(lambda x: decay(x,*poptA1)-2, 16)

# error propagation for 2=A/(x0**tau)
sigma_A, sigma_tau = np.sqrt(np.diag(pcovA1))
A, tau = poptA1
varx0 = (np.log(A/2)/(tau*np.log(tau)**2))**2 * sigma_tau**2 + \
(1/(A*np.log(tau)))**2 * sigma_A**2
sigmax0 = np.sqrt(varx0)

display(Latex("We find that $T_c={:.0f} \pm {:.0f}.$".format(x0,sigmax0)))


We find that $T_c=11 \pm 133.$

It's difficult to estimate the correctness of our model when our uncertainties are this large, but even if we found a better model, the uncertainties are simply too large to meaningfully find $T_c$ using this method. We get no information from this investigation, and we use our previous conclusion that $T_c$ lies somewhere between $10K$ and $20K$.