Companion to 'Coding closed and open quantum systems in MATLAB' (arXiv:1911.04906)
Runnable numbl scripts for the article's examples (Norambuena, Tancara, Coto),
reproduced verbatim, plus a README guide and GitHub Pages deploy workflow.
15 changed files+653−0
.github/workflows/deploy.ymladded+44−0View file
@@ -0,0 +1,44 @@
1+name: Deploy numbl project to GitHub Pages
2+
3+# Bundles this numbl project with the browser IDE and publishes it to GitHub
4+# Pages on every push to main.
5+#
6+# One-time setup: Settings → Pages → "Build and deployment" → Source →
7+# "GitHub Actions".
8+
9+on:
10+ push:
11+ branches: [main]
12+ workflow_dispatch:
13+
14+permissions:
15+ contents: read
16+ pages: write
17+ id-token: write
18+
19+concurrency:
20+ group: pages
21+ cancel-in-progress: false
22+
23+jobs:
24+ build:
25+ runs-on: ubuntu-latest
26+ steps:
27+ - uses: actions/checkout@v4
28+ - uses: flatironinstitute/numbl/.github/actions/build-site@main
29+ with:
30+ project-dir: .
31+ # Build numbl from the main branch (development version) rather than
32+ # the published npm release — these examples need recent fixes
33+ # (complex eig, trace, parfor, box, legend).
34+ numbl-ref: main
35+
36+ deploy:
37+ needs: build
38+ runs-on: ubuntu-latest
39+ environment:
40+ name: github-pages
41+ url: ${{ steps.deployment.outputs.page_url }}
42+ steps:
43+ - id: deployment
44+ uses: actions/deploy-pages@v4
.gitignoreadded+4−0View file
@@ -0,0 +1,4 @@
1+# Local output from `numbl build-site`
2+/_site/
3+/dist/
4+node_modules/
CavityQEDTransitionPhase.madded+72−0View file
@@ -0,0 +1,72 @@
1+Nsim = 25; % Number of simulations
2+L = 2; % Number of cavities
3+wc = 1; % Cavity frequency
4+g = 1e-2*wc; % Atom-light coupling
5+J = 1e-4*wc; % Coupling between cavities
6+Nph = 2; % Number of photons per cavity
7+dimFock = Nph+1; % Dimension Fock space for photons
8+dimT = 2*dimFock; % Dimension atom+cavity system
9+Deltai = 10^(-2)*g; % Initial detuning
10+Deltaf = 10^(+2)*g; % Final detuning
11+xi = log10(Deltai/g); % Initial detuning in log scale
12+xf = log10(Deltaf/g); % Final detuning in log scale
13+dx = (xf-xi)/(Nsim-1); % Step dx
14+x = xi:dx:xf; % Vector x to plot transition phase
15+OP_JC = zeros(size(x)); % Order parameter Jaynes-Cummings model
16+OP_R = zeros(size(x)); % Order parameter Rabi model
17+A = cell(1,L); % Cell array to storage \hat{a}_i operators
18+Sp = cell(1,L); % Cell array to storage \hat{\sigma}_i^+ operators
19+N_ex = cell(1,L); % Cell array to storage N_i operators
20+Iatom = eye(2); % Identity matrix atom system
21+Icav = eye(dimFock ); % Identity matrix cavity system
22+Is = eye(2*dimFock); % Identity matrix atom+cavity system
23+for i=1:L
24+ A{i} = acav(i,L,Nph,Is,Iatom); % Photon operator \hat{a}_i
25+ Sp{i} = sigmap(i,L,Is,Icav); % Atom operator \hat{\sigma}_i
26+ N_ex{i} = A{i}'*A{i}+Sp{i}*Sp{i}'; % number of excitations per cavity
27+end
28+Ad = ones(L); % Adyacent matrix
29+Ad = triu(Ad)-eye(L); % A(i,j)=1 for j>i
30+Hhopp = zeros(dimT^L,dimT^L); % Interaction Hamiltonian
31+for i=1:L
32+ for j=1:L
33+ Hhopp = Hhopp - J*Ad(i,j)*A{i}'*A{j} - J*Ad(i,j)*A{i}*A{j}';
34+ end
35+end
36+Nt = 10000; % Length time vector
37+ti = 0.01/J; % Initial time
38+tf = 1/J; % Final time
39+dt = (tf-ti)/(Nt-1); % Step time dt
40+t = ti:dt:tf; % Time vector
41+Lambda_R = zeros(Nsim,Nt);
42+Lambda_JC = zeros(Nsim,Nt);
43+parfor n=1:Nsim
44+ D = g*10^(x(n)); % Detuning in each simulation
45+ Model = 'Rabi';
46+ [OP_R(n), P1m_R] = QuantumSimulationCavityArray(wc,D,g,Nph,L,A,Sp,N_ex,Hhopp,t,Model);
47+ Model = 'Jaynes-Cummings';
48+ [OP_JC(n),P1m_JC] = QuantumSimulationCavityArray(wc,D,g,Nph,L,A,Sp,N_ex,Hhopp,t,Model);
49+ Lambda_R(n,:) = -1/L*log2(P1m_R); % Rate function for the Rabi model
50+ Lambda_JC(n,:) = -1/L*log2(P1m_JC); % Rate function for the Jaynes-Cummings model
51+end
52+figure()
53+box on
54+hold on
55+plot(J*t,Lambda_R(end,:),'b-','Linewidth',3)
56+plot(J*t,Lambda_JC(end,:),'r-','Linewidth',3)
57+hold off
58+xlabel('$Jt$','Interpreter','LaTex','Fontsize', 30)
59+ylabel('$\Lambda(t)$','Interpreter','LaTex','Fontsize', 30)
60+set(gca,'fontsize',21)
61+legend({'$\mbox{RH}$','$\mbox{JCH}$'},'Interpreter','latex','Fontsize', 21,'Location','best')
62+
63+figure()
64+box on
65+hold on
66+plot(x,real(OP_R),'.b','Markersize',30)
67+plot(x,real(OP_JC),'.r','Markersize',30)
68+hold off
69+xlabel('$\mbox{Log}_{10}(\Delta/g)$','Interpreter','LaTex','Fontsize', 30)
70+ylabel('$\mbox{OP}$','Interpreter','LaTex','Fontsize', 30)
71+set(gca,'fontsize',21)
72+legend({'$\mbox{RH}$','$\mbox{JCH}$'},'Interpreter','latex','Fontsize', 21,'Location','best')
Example1.madded+29−0View file
@@ -0,0 +1,29 @@
1+J = 1; % Value of the coupling J
2+B = 0.1*J; % Value of the magnetic field
3+Sx = [0 1;1 0]; % S_x operator
4+Sz = [1 0; 0 -1]; % S_z operator
5+I = eye(2); % Identity matrix
6+Hspins = -J*kron(Sx,Sx)-B*(kron(Sx,I)+kron(I,Sx)); % Hamiltonian
7+down = [0 1]'; % Quantum state down = [0 1]^T
8+Psi_0 = kron(down,down); % Initial wavefunction
9+T = 2*pi/(2*B); % Period of time
10+Nt = 1000; % Number of steps to construct time vector
11+ti = 0; % Initial time
12+tf = 2*T; % Final time
13+dt = (tf-ti)/(Nt-1); % Step time dt
14+t = ti:dt:tf; % Time vector
15+U = expm(-1i*Hspins*dt); % Time propagator operator U(dt)
16+SSz = (kron(Sz,I)+kron(I,Sz))/2; % Operator S1^z + S2^z
17+Mz = zeros(size(t)); % Average magnetization
18+for n=1:length(t) % Iteration to find Psi(t) and Mz(t)
19+ if n==1
20+ Psi = Psi_0; % Initial wavefunction
21+ else
22+ Psi = U*Psi; % Wavefuntion at time t_n
23+ end
24+ Mz(n) = Psi'*SSz*Psi; % Average magnetization at time t_n
25+end
26+plot(t/T,real(Mz),'r-','LineWidth',3)
27+xlabel('$t/T$','Interpreter','LaTex','Fontsize', 30)
28+ylabel('$\langle M_z \rangle$','Interpreter','LaTex','Fontsize', 30)
29+set(gca,'fontsize',21)
Example2.madded+51−0View file
@@ -0,0 +1,51 @@
1+dim = 4; % Dimension of the total Hilbert space
2+Is = eye(dim); % Identity matrix for the total Hilbert space
3+J = 1; % Value of the coupling term J
4+B = 0.1*J; % Value of the magnetic field B
5+Sx = [0 1;1 0]; % S_x operator for one spin
6+Sz = [1 0;0 -1]; % S_z operator for one spin
7+I = eye(2); % Identity matrix for one spin 1/2
8+H = -J*kron(Sx,Sx)-B*(kron(Sx,I)+kron(I,Sx)); % Hamiltonian of the system
9+L_H = -1i*kron(Is,H)+1i*kron(H.',Is); % Lindblad operator L_H
10+gamma_1 = 0.1*B; % Decay rate gamma_1
11+gamma_2 = 0.5*B; % Decay rate gamma_2
12+S_minus =[0 0; 1 0]; % Lowering operator of the particle 1
13+L1 = kron(S_minus,I); % Lowering operator of the particle 1 in the total Hilbert space
14+DL_1 = gamma_1*(kron(conj(L1),L1)-0.5*kron(Is,L1'*L1)-0.5*kron(L1.'*conj(L1),Is));
15+L2 = kron(I,S_minus); % Lowering operator of the particle 2 in the total Hilbert space
16+DL_2 = gamma_2*(kron(conj(L2),L2)-0.5*kron(Is,L2'*L2)-0.5*kron(L2.'*conj(L2),Is));
17+L_diss = DL_1 + DL_2; % Lindbladian L_diss
18+L = L_H + L_diss; % Total Lindblad operator
19+TOL = 1e-10;
20+[R_sort,L_sort,lambda_sort] = sortingEigenvalues(dim,TOL,L);
21+
22+down = [0 1]'; % Quantum state down = [0 1]
23+Psi_0 = kron(down,down); % Initial wavefunction
24+rho_0 = Psi_0*Psi_0'; % Initial density matrix
25+Nt = 3000; % Number of steps for time
26+T = 2*pi/(2*B) ; % Period of time
27+ti = 0; % Intial time
28+tf = 4*T; % Final time
29+dt = (tf-ti)/(Nt-1); % Step for time
30+t = ti:dt:tf; % Time vector
31+Mz = zeros(size(t)); % Initial average magnetization
32+SSz = (kron(Sz,I)+kron(I,Sz))/2; % Operator S1^z+S2^z
33+for n=1:length(t) % Iteration to find general sulution of rho(t)
34+ rho = zeros(dim,dim);
35+ for k=1:length(lambda_sort)
36+ Lk = L_sort{k};
37+ Rk = R_sort{k};
38+ ck = trace(rho_0*Lk);
39+ rho = rho + ck*exp(lambda_sort(k)*t(n))*Rk; % General sulution for rho(t)
40+ end
41+ Mz(n) = trace(SSz*rho); % General sulution for Mz(t)
42+end
43+figure()
44+hold on
45+plot(t/T,real(Mz),'r-','LineWidth',3)
46+plot(t/T, exp(real(lambda_sort(end))*t),'k--','LineWidth',2)
47+plot(t/T,-exp(real(lambda_sort(end))*t),'k--','LineWidth',2)
48+hold off
49+xlabel('$Bt$','Interpreter','LaTex','Fontsize', 30)
50+ylabel('$\langle M_z \rangle$','Interpreter','LaTex','Fontsize', 30)
51+set(gca,'fontsize',21)
IsingModelTransitionPhase_new.madded+75−0View file
@@ -0,0 +1,75 @@
1+Nspins = 6; % Number of spins
2+H0 = zeros(2^Nspins,2^Nspins);
3+H1 = zeros(2^Nspins,2^Nspins);
4+J = 0;
5+alpha = 0.2; % Parameter \alpha
6+for j=1:Nspins
7+ for i=j+1:Nspins
8+ Jij = abs(i-j)^(-alpha);
9+ J = J + (Nspins-1)^(-1)*Jij;
10+ end
11+end
12+B = J/0.42; % Magnetic field
13+Sx = [0 1;1 0];
14+Sz = [1 0;0 -1];
15+for i=1:Nspins
16+ Szi = getSci(Sz,i,Nspins);
17+ H0 = H0 - B*Szi; % Hamiltonian for the magnetic field
18+ for j=1:Nspins
19+ if i~=j
20+ Sxi = getSci(Sx,i,Nspins);
21+ Sxj = getSci(Sx,j,Nspins);
22+ Vij = abs(i-j)^(-alpha)/J;
23+ H1 = H1 - Vij*Sxi*Sxj; % Interaction Hamiltonian
24+ end
25+ end
26+end
27+H = H0 + H1; % Total hamiltonian
28+xr = [1 1]'/sqrt(2); % Single-particle state |Psi_{-->}>
29+xl = [-1 1]'/sqrt(2); % Single-particle state |Psi_{<--}>
30+Xr = xr;
31+Xl = xl;
32+for n=1:Nspins-1
33+ Xr = kron(Xr,xr); % Many-body state |Psi_{-->}>
34+ Xl = kron(Xl,xl); % Many-body state |Psi_{<--}>
35+end
36+PSI_0 = Xr; % Initial condition
37+ti = 0; % Initial time
38+tf = 22; % Final time
39+Nt = 10000; % Number of steps
40+dt = (tf-ti)/(Nt-1); % Step time dt
41+t = ti:dt:tf; % Time vector
42+U = expm(-1i*H*dt); % Time propagator operator U(dt)
43+Mx = zeros(size(t)); % Average Magnetization <M_x(t)>
44+Lambda = zeros(size(t)); % Rate function \Lambda(t)
45+SSx = 0;
46+for i=1:Nspins
47+ Sxi = getSci(Sx,i,Nspins);
48+ SSx = SSx + Sxi/Nspins; % Magnetization operator
49+end
50+for n=1:length(t)
51+ if n==1
52+ PSI = PSI_0; % Initial wavefunction
53+ else
54+ PSI = U*PSI; % Wavefunction at time t_n
55+ end
56+ Pr = abs(Xr'*PSI)^2; % Probability state |Psi_{-->}>
57+ Pl = abs(Xl'*PSI)^2; % Probability state |Psi_{<--}>
58+ Lambda(n) = min(-Nspins^(-1)*log(Pr),-Nspins^(-1)*log(Pl));
59+ Mx(n) = PSI'*SSx*PSI; % Average magnetization along x-axis
60+end
61+
62+figure()
63+plot(B*t,Lambda,'b-','Linewidth',3)
64+xlabel('$B t$','Interpreter','LaTex','Fontsize', 30)
65+ylabel('$\Lambda(t)$','Interpreter','LaTex','Fontsize', 30)
66+set(gca,'fontsize',21)
67+xlim([0 5])
68+
69+figure()
70+box on
71+plot(t*B,real(Mx),'r-','Linewidth',2)
72+xlabel('$B t$','Interpreter','LaTex','Fontsize', 30)
73+ylabel('$\langle M_x\rangle$','Interpreter','LaTex','Fontsize', 30)
74+set(gca,'fontsize',21)
75+xlim([0 100])
NonMarkovianDynamicsPureDephasing.madded+106−0View file
@@ -0,0 +1,106 @@
1+sz = [1 0; % Pauli matrix s_z
2+ 0 -1];
3+Nw = 5000; % Number of points for \omega
4+wi = 0.01; % Initial frequency
5+wf = 5; % Final frequency
6+dw = (wf-wi)/(Nw-1); % Step frequency d\omega
7+w = wi:dw:wf; % Frequency vector
8+alpha = 5; % Strength spectral density function J1(w)
9+s = 2.5; % Ohmic-parameter s
10+wc = 0.1; % Cut-off frequency
11+J1 = alpha*wc^(1-s)*w.^s.*exp(-w/wc); % Spectral density funcion J1(w)
12+J0 = 0.2; % Strength spectral density function J2(w)
13+w0 = 2; % Resonant frequency
14+Gamma = 0.1; % Width of the spectral density function J2(w)
15+J2 = J0*(Gamma/2)./((w-w0).^2+(Gamma/2)^2).*w.^s./(w/w0+1).^2; % Spectral density funcion J2(w)
16+figure
17+box on
18+hold on
19+plot(w,J1,'r-','Linewidth',2)
20+plot(w,J2,'b-','Linewidth',2)
21+hold off
22+xlabel('$\omega$','Interpreter','LaTex','Fontsize', 30)
23+ylabel('$J(\omega)$','Interpreter','LaTex','Fontsize', 30)
24+legend({'$J_1(\omega)$','$J_2(\omega)$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
25+set(gca,'fontsize',21)
26+xlim([0 5])
27+Nt = 5000; % Number of points for time
28+ti = 0; % Initial time
29+tf = 100; % Final time
30+dt = (tf-ti)/(Nt-1); % Step time dt
31+t = ti:dt:tf; t= t'; % Time vector
32+Psi_0 = [1 1]'/sqrt(2); % Initial wavefunction
33+rho1 = Psi_0*Psi_0'; % Initial density matrix
34+rho2 = rho1;
35+p11 = zeros(size(t));
36+p22 = zeros(size(t));
37+p12 = zeros(size(t));
38+p21 = zeros(size(t));
39+wa = ones(size(t))*w; % Auxiliar frecuency vector
40+J1a = ones(size(t))*J1; % Auxiliar J1 vector
41+J2a = ones(size(t))*J2; % Auxiliar J2 vector
42+ta = t*ones(size(w)); % Auxiliar time vector
43+T = 0.001*w0; % Temperature
44+gamma1 = sum(J1a./wa.*sin(wa.*ta).*coth(wa/T/2),2)*dw; % Rate gamma_1(t)
45+gamma_1_teo = alpha*wc*gamma(s)*sin(s*atan(wc*t))./(1+(wc*t).^2).^(s/2);
46+gamma2 = sum(J2a./wa.*sin(wa.*ta).*coth(wa/T/2),2)*dw; % Rate gamma_2(t)
47+figure()
48+box on
49+hold on
50+plot(t,gamma1,'r-','Linewidth',2)
51+plot(t,gamma_1_teo,'k--','Linewidth',2)
52+plot(t,gamma2,'b-','Linewidth',2)
53+hold off
54+xlabel('$t$','Interpreter','LaTex','Fontsize', 30)
55+ylabel('$\gamma(t)$','Interpreter','LaTex','Fontsize', 30)
56+legend({'$\gamma_1(t)$','$\gamma_1^{\rm teo}(t)$','$\gamma_2(t)$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
57+set(gca,'fontsize',21)
58+xlim([0 100])
59+nk = 15; % nk steps per interval dt
60+C1 = zeros(size(t));
61+C2 = zeros(size(t));
62+N1 = zeros(size(t));
63+N2 = zeros(size(t));
64+for n = 1:length(t)
65+ for k = 1:nk
66+ L1 = gamma1(n)*(sz*rho1*sz -rho1); % Lindbladian for gamma_1
67+ L2 = gamma2(n)*(sz*rho2*sz -rho2); % Lindbladian for gamma_2
68+ rho1_pred = rho1 + L1*dt/nk; % Predictor \rho_1
69+ rho2_pred = rho2 + L2*dt/nk; % Predictor \rho_2
70+ rho1_m = 0.5*(rho1+rho1_pred); % Mean \rho_1
71+ rho2_m = 0.5*(rho1+rho2_pred); % Mean \rho_1
72+ L1 = gamma1(n)*(sz*rho1_m*sz -rho1_m); % Lindbladian using mean \rho_1
73+ L2 = gamma2(n)*(sz*rho2_m*sz -rho2_m); % Lindbladian using mean \rho_2
74+ rho1 = rho1 + L1*dt/nk; % Density matrix \rho_1
75+ rho2 = rho2 + L2*dt/nk; % Density matrix \rho_2
76+ end
77+ C1(n) = 2*abs(rho1(1,2)); % Coherence C1(t)
78+ C2(n) = 2*abs(rho2(1,2)); % Coherence C2(t)
79+ f1 = (abs(gamma1(1:n))-gamma1(1:n));
80+ N1(n) = sum(f1)*dt; % Degree of non-Markovianity for gamma_1
81+ f2 = (abs(gamma2(1:n))-gamma2(1:n));
82+ N2(n) = sum(f2)*dt; % Degree of non-Markovianity for gamma_2
83+end
84+figure()
85+box on
86+hold on
87+plot(t,C1,'r-','Linewidth',2)
88+plot(t,C2,'b-','Linewidth',2)
89+hold off
90+xlabel('$t$','Interpreter','LaTex','Fontsize', 30)
91+ylabel('$C(t)$','Interpreter','LaTex','Fontsize', 30)
92+legend({'$C_1(t)$','$C_2(\omega)$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
93+set(gca,'fontsize',21)
94+xlim([0 50])
95+
96+figure()
97+box on
98+hold on
99+plot(t,N1,'r-','Linewidth',2)
100+plot(t,N2,'b-','Linewidth',2)
101+hold off
102+xlabel('$t$','Interpreter','LaTex','Fontsize', 30)
103+ylabel('$\mathcal{N}_{\gamma}(t)$','Interpreter','LaTex','Fontsize', 30)
104+legend({'$\mathcal{N}_{\gamma_1}(t)$','$\mathcal{N}_{\gamma_2}(t)$'},'Interpreter','latex','Fontsize', 21,'Location','best')
105+set(gca,'fontsize',21)
106+xlim([0 100])
QuantumSimulationCavityArray.madded+40−0View file
@@ -0,0 +1,40 @@
1+function [OP,P1m]= QuantumSimulationCavityArray(wc,D,g,Nphoton,L,A,Sp,N_ex,Hhopp,t,Model)
2+HJC = zeros((Nphoton+1)^L*2^L,(Nphoton+1)^L*2^L);
3+for i=1:L
4+ switch Model
5+ case 'Jaynes-Cummings'
6+ HJC = HJC + wc*A{i}'*A{i} + (D+wc)*Sp{i}*Sp{i}' + g*(Sp{i}*A{i}+Sp{i}'*A{i}');
7+ case 'Rabi'
8+ HJC = HJC + wc*A{i}'*A{i} + (D+wc)*Sp{i}*Sp{i}' + g*(Sp{i}+Sp{i}')*(A{i}+A{i}');
9+ end
10+end
11+H = HJC + Hhopp; % Total Hamiltonian
12+dt = t(2)-t(1); % Step time dt
13+U = expm(-1i*H*dt); % Time propagator with step dt
14+up = [1 0]'; % Excited state for the two-level system
15+down = [0 1]'; % Ground state for the two-level system
16+Fock = eye(Nphoton+1); % Fock states
17+theta1 = 0.5*atan(2*g*sqrt(1)/D);
18+phi_1m = cos(theta1)*kron(down,Fock(:,2))...
19+ -sin(theta1)*kron(up,Fock(:,1)); % Initial state |1,->
20+PSI_0 = phi_1m;
21+for k=1:L-1
22+ PSI_0 = kron(PSI_0,phi_1m); % Many body wavefunction |1,->...|1,-> (L terms)
23+end
24+dn_T = zeros(size(t)); % Standard deviation dn_i = <n_i^2>-<n_i>^2
25+P1m = zeros(size(t)); % Population |1,-> at time t
26+for n=1:length(t)
27+ if n==1
28+ PSI = PSI_0; % Initial wavefunction
29+ else
30+ PSI = U*PSI; % Wavefunction at time t_n
31+ end
32+ dn = 0;
33+ for i=1:L
34+ dn = dn + PSI'*N_ex{i}^2*PSI-(PSI'*N_ex{i}*PSI)^2;
35+ end
36+ P1m(n) = abs(PSI_0'*PSI)^2; % Ground state probability
37+ dn_T(n) = dn; % Standard deviation of the number of excitations
38+end
39+OP = mean(dn_T); % Order parameter
40+end
README.mdadded+78−0View file
@@ -0,0 +1,78 @@
1+# Companion to Closed & Open Quantum Systems
2+
3+A companion to the tutorial article **"Coding closed and open quantum systems in
4+MATLAB: applications in quantum optics and condensed matter"** by Ariel
5+Norambuena, Diego Tancara, and Raúl Coto —
6+[arXiv:1911.04906](https://arxiv.org/abs/1911.04906).
7+
8+The article presents MATLAB code for simulating many-body quantum systems:
9+magnetization dynamics of the closed and open Ising model, a dynamical quantum
10+phase transition in cavity-QED arrays, Markovian dynamics of interacting
11+two-level systems, and the non-Markovian pure-dephasing spin-boson model. This
12+project collects those programs as runnable [numbl](https://numbl.org) scripts
13+so you can execute and edit them in your browser — no MATLAB, no install. Pick a
14+script in the panel on the left and click **Run**. Edits stay in memory for your
15+session.
16+
17+> This README is a guide to *which* program to open, not a replacement for the
18+> article. Read the article for the physics, the equations, and the derivations;
19+> come here to run the code. The scripts are reproduced **verbatim** from the
20+> article's source (arXiv ancillary files) and run unmodified. The only
21+> repackaging: the article presents Example 1 and Example 2 as a sequence of
22+> code listings, which are concatenated here into single runnable
23+> [`Example1.m`](Example1.m) and [`Example2.m`](Example2.m) files.
24+
25+## Closed quantum systems (article §2)
26+
27+| Script | Article | What it computes |
28+| --- | --- | --- |
29+| [`Example1.m`](Example1.m) | §2.1 | Magnetization dynamics of the closed two-spin Ising model via the time propagator `expm(-iHt)` |
30+| [`IsingModelTransitionPhase_new.m`](IsingModelTransitionPhase_new.m) | §2.1 | Dynamical quantum phase transition of the Ising model (6 spins): rate function and magnetization |
31+| [`CavityQEDTransitionPhase.m`](CavityQEDTransitionPhase.m) | §2.2 | Dynamical phase transition in a cavity-QED array — Jaynes–Cummings vs. Rabi (order parameter and rate function) |
32+
33+## Open quantum dynamics (article §3)
34+
35+| Script | Article | What it computes |
36+| --- | --- | --- |
37+| [`Example2.m`](Example2.m) | §3.1 | Markovian master equation for two interacting qubits (open Ising), solved by diagonalizing the Lindbladian |
38+| [`TwoLevelSystemCoupledLight.m`](TwoLevelSystemCoupledLight.m) | §3.2 | A two-level atom coupled to a photon reservoir; numerical solution compared to the exact one |
39+| [`NonMarkovianDynamicsPureDephasing.m`](NonMarkovianDynamicsPureDephasing.m) | §3.4 | Non-Markovian pure-dephasing spin-boson model: time-dependent rates, coherence, and a non-Markovianity measure |
40+
41+## Helper functions
42+
43+These are called by the scripts above; they are functions, not standalone
44+programs, so there is nothing to "Run" directly.
45+
46+| File | Used by | Builds |
47+| --- | --- | --- |
48+| [`getSci.m`](getSci.m) | `IsingModelTransitionPhase_new.m` | the many-body Pauli operator on site *i* |
49+| [`acav.m`](acav.m) | `CavityQEDTransitionPhase.m` | the cavity annihilation operator on site *i* |
50+| [`sigmap.m`](sigmap.m) | `CavityQEDTransitionPhase.m` | the atomic raising operator on site *i* |
51+| [`QuantumSimulationCavityArray.m`](QuantumSimulationCavityArray.m) | `CavityQEDTransitionPhase.m` | the Jaynes–Cummings/Rabi Hamiltonian, time evolution, and observables |
52+| [`sortingEigenvalues.m`](sortingEigenvalues.m) | `Example2.m`, `TwoLevelSystemCoupledLight.m` | the sorted left/right eigenmatrices of the Lindbladian |
53+
54+## Notes
55+
56+- **Time evolution** of closed systems uses the matrix exponential
57+ `U = expm(-1i*H*dt)`. **Open systems** are solved by diagonalizing the
58+ Lindblad superoperator with `eig` and reconstructing the density matrix from
59+ its eigenmodes (`sortingEigenvalues`).
60+- `CavityQEDTransitionPhase.m` is the heaviest example — it sweeps 25 detunings
61+ for two models with long time evolutions, so it takes a couple of minutes in
62+ the browser. The others run in a few seconds.
63+- The article uses `parfor` in the cavity-QED sweep; numbl runs it as an
64+ ordinary (serial) loop, which gives the same result.
65+
66+## How it's deployed
67+
68+On every push to `main`, the workflow in
69+[`.github/workflows/deploy.yml`](.github/workflows/deploy.yml) bundles these
70+files with the numbl browser IDE and publishes the result to GitHub Pages — no
71+server, nothing to build by hand. Edit `numbl-project.json` to change the title
72+or which file opens first.
73+
74+## Reference
75+
76+Ariel Norambuena, Diego Tancara, Raúl Coto, *Coding closed and open quantum
77+systems in MATLAB: applications in quantum optics and condensed matter*,
78+[arXiv:1911.04906](https://arxiv.org/abs/1911.04906).
TwoLevelSystemCoupledLight.madded+69−0View file
@@ -0,0 +1,69 @@
1+Omega = 1; % Raby frequency from |1> to |2>
2+gamma0 = 0.2*Omega; % Decay rate
3+dim = 2; % Dimension Hilbert space two-level system
4+Is = eye(dim); % Identity matrix Hilbert space
5+v1 = Is(:,1); % Excited state for the atom
6+v2 = Is(:,2); % Ground state fro the atom
7+s11 = v1*v1'; % Atom operator sigma_{11}
8+s22 = v2*v2'; % Atom operator sigma_{22}
9+sp = v1*v2'; % Atom operator sigma_{+}
10+sm = v2*v1'; % Atom operator sigma_{-}
11+HL = -0.5*Omega*(sp+sm); % Hamiltonian of the two-level system
12+N = 0; % Mean number of photons at zero temperature
13+Lrad = gamma0*(N+1)*(kron(conj(sm),sm)-0.5*kron(Is,sm'*sm)-0.5*kron(sm.'*conj(sm),Is)) + ...
14+ gamma0*N*(kron(conj(sp),sp)-0.5*kron(Is,sp'*sp)-0.5*kron(sp.'*conj(sp),Is));
15+L = -1i*kron(Is,HL)+1i*kron(HL.',Is)+Lrad; % Lindblad operator
16+TOL = 1e-6;
17+[R_sort,L_sort,lambda_sort] = sortingEigenvalues(dim,TOL,L);
18+psi_0 = v2; % Ground state
19+rho_0 = psi_0*psi_0'; % Initial density matrix
20+Nt = 100000; % Number of steps for time
21+ti = 0; % Initial time
22+tf = 50/Omega; % Final time
23+dt = (tf-ti)/(Nt-1); % Step time dt
24+t = ti:dt:tf; % Time vector
25+rho11 = zeros(size(t)); % Matrix elements \rho_{ij} = <i|\rho|j>
26+sigmap = zeros(size(t));
27+for n=1:Nt % General solution
28+ rho = zeros(dim,dim);
29+ for k=1:length(lambda_sort)
30+ Lk = L_sort{k};
31+ Rk = R_sort{k};
32+ ck = trace(rho_0*Lk);
33+ rho = rho + ck*exp(lambda_sort(k)*t(n))*Rk;
34+ end
35+ rho11(n) = rho(1,1);
36+ sigmap(n) = trace(rho*sp);
37+end
38+% Exact Solution at zero temperature
39+% p_e: population excited state, sp = <\sigma_x>
40+mu = 1i*sqrt((gamma0/4)^2-Omega^2);
41+pe_exact = Omega^2/(gamma0^2+2*Omega^2)*(1-exp(-3*gamma0*t/4).*(cos(mu*t)+3*gamma0/4/mu*sin(mu*t)));
42+if gamma0~=0
43+ sp_exact = -1i*Omega*gamma0/(gamma0^2+2*Omega^2)*(1-exp(-3*gamma0*t/4).*(cos(mu*t)+(gamma0/4/mu)-Omega^2/gamma0/mu*sin(mu*t)));
44+else
45+ sp_exact = -1i/2*sin(Omega*t);
46+end
47+figure
48+box on
49+hold on
50+plot(t*Omega,real(rho11),'r-','Linewidth',2)
51+plot(t*Omega,pe_exact,'b--','Linewidth',2)
52+hold off
53+xlabel('$\Omega t$','Interpreter','LaTex','Fontsize', 30)
54+ylabel('$p_e(t)$','Interpreter','LaTex','Fontsize', 30)
55+legend({'$\mbox{Numerical}$','$\mbox{Exact}$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
56+set(gca,'fontsize',21)
57+xlim([0 50])
58+
59+figure
60+box on
61+hold on
62+plot(t*Omega,imag(sigmap),'r-','Linewidth',2)
63+plot(t*Omega,imag(sp_exact),'b--','Linewidth',2)
64+hold off
65+xlabel('$\Omega t$','Interpreter','LaTex','Fontsize', 30)
66+ylabel('$\mbox{Im}\langle \hat{\sigma}_+ \rangle$','Interpreter','LaTex','Fontsize', 30)
67+legend({'$\mbox{Numerical}$','$\mbox{Exact}$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
68+set(gca,'fontsize',21)
69+xlim([0 50])
acav.madded+12−0View file
@@ -0,0 +1,12 @@
1+function x = acav(i,L,Nphoton,Is,Iatom)
2+a = diag(sqrt(1:Nphoton)',1);
3+a = kron(Iatom,a);
4+Op_total = cell(1,L);
5+for site = 1:L
6+ Op_total{site} = Is+double(eq(i,site))*(a-Is);
7+end
8+x = Op_total{1};
9+for site = 2:L
10+ x = kron(x,Op_total{site});
11+end
12+end
getSci.madded+11−0View file
@@ -0,0 +1,11 @@
1+function Sci = getSci(sc,i,Nspins)
2+Is = eye(2);
3+Op_total = cell(1,Nspins);
4+for site = 1:Nspins
5+ Op_total{site} = Is+double(eq(i,site))*(sc-Is);
6+end
7+Sci = Op_total{1};
8+for site = 2:Nspins
9+ Sci = kron(Sci,Op_total{site});
10+end
11+end
numbl-project.jsonadded+4−0View file
@@ -0,0 +1,4 @@
1+{
2+ "title": "Companion to Closed & Open Quantum Systems",
3+ "entry": "README.md"
4+}
sigmap.madded+14−0View file
@@ -0,0 +1,14 @@
1+function x = sigmap(i,L,Is,Icav)
2+up = [1 0]';
3+down = [0 1]';
4+sigma_p = up*down';
5+sigma_p = kron(sigma_p,Icav);
6+Op_total = cell(1,L);
7+for site = 1:L
8+ Op_total{site} = Is+double(eq(i,site))*(sigma_p-Is);
9+end
10+x = Op_total{1};
11+for site = 2:L
12+ x = kron(x,Op_total{site});
13+end
14+end
sortingEigenvalues.madded+44−0View file
@@ -0,0 +1,44 @@
1+function [R_sort,L_sort,lambda_sort] = sortingEigenvalues(dim,TOL,L)
2+
3+[R,DR] = eig(L); % Right eigenvectors and eigenvalues
4+[L,DL] = eig(L'); % Left eigenvectors and eigenvalues
5+eig_R = diag(DR); % Right eigenvalues written as a vector
6+eig_L = diag(DL); % Left eigenvalues written as a vector
7+ind_RL = zeros(dim*dim,2);
8+count = 1;
9+for n=1:dim*dim % Sorting of eigenvalues
10+ an = eig_R(n);
11+ for m=1:dim*dim
12+ bm = eig_L(m);
13+ if(abs(real(an)-real(bm))<TOL && abs(imag(an)-imag(bm))<TOL && count<=dim*dim)
14+ ind_RL(count,1) = n;
15+ ind_RL(count,2) = m;
16+ count = count + 1;
17+ end
18+ end
19+end
20+eig_L = eig_L(ind_RL(:,2)'); % Final sorting
21+eig_R = eig_R(ind_RL(:,1)');
22+L = L(:,ind_RL(:,2)');
23+R = R(:,ind_RL(:,1)');
24+lambda = eig_R;
25+[~,ind] = sort(lambda); % Sorting of eigenvalues
26+lambda_sort = lambda(ind); % \lambda_k eigenvalues
27+L = L(:,ind);
28+R = R(:,ind);
29+
30+% Final R_sort and L_sort matrices
31+R_sort = cell(1,length(lambda_sort));
32+L_sort = cell(1,length(lambda_sort));
33+for k=1:length(lambda_sort)
34+ R_sort{k} = reshape(R(:,k),dim,dim);
35+ L_sort{k} = reshape(L(:,k),dim,dim);
36+ Rk = R_sort{k};
37+ Lk = L_sort{k};
38+ Ck = trace(Lk*Rk);
39+ Lk = Lk/sqrt(Ck); % Normalized left eigenmatrices
40+ Rk = Rk/sqrt(Ck); % Normalized right eigenmatrices
41+ R_sort{k} = Rk;
42+ L_sort{k} = Lk;
43+end
44+end