Quantum optics in numbl: companion to arXiv:2309.14354
Runnable numbl scripts for each example in Nilakantha Meher's tutorial
'Quantum optics in MATLAB' (arXiv:2309.14354), one .m file per code block,
plus a README guide and GitHub Pages deploy workflow.
21 changed files+501−0
.github/workflows/deploy.ymladded+40−0View file
@@ -0,0 +1,40 @@
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+
32+ deploy:
33+ needs: build
34+ runs-on: ubuntu-latest
35+ environment:
36+ name: github-pages
37+ url: ${{ steps.deployment.outputs.page_url }}
38+ steps:
39+ - id: deployment
40+ 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/
README.mdadded+84−0View file
@@ -0,0 +1,84 @@
1+# Quantum Optics in numbl
2+
3+A companion to the tutorial article **"Quantum optics in MATLAB"** by Nilakantha
4+Meher — [arXiv:2309.14354](https://arxiv.org/abs/2309.14354).
5+
6+The article teaches the basics of theoretical quantum optics through short,
7+self-contained MATLAB scripts: building quantum states and operators as
8+matrices, then using them to compute physical properties and time dynamics.
9+This project collects those examples as runnable
10+[numbl](https://numbl.org) scripts so you can execute and edit them directly in
11+your browser — no MATLAB, no install. Pick a script in the panel on the left and
12+click **Run**. Edits stay in memory for your session.
13+
14+> This README is a guide to *which* example to open, not a replacement for the
15+> article. Read the article for the physics, the equations, and the discussion;
16+> come here to run the code. Each numbered script corresponds to one code block
17+> in the article, in the same order, and is annotated with its section number.
18+> The code is transcribed faithfully from the article; the only additions are
19+> axis labels and legends on the plots.
20+
21+## Quantum states (Sec. II)
22+
23+| Script | Article | What it shows |
24+| --- | --- | --- |
25+| [`ex01_number_states.m`](ex01_number_states.m) | II.A | Number states `\|0>, \|1>, \|2>` as columns of the identity matrix |
26+| [`ex02_number_state_20.m`](ex02_number_state_20.m) | II.A | A higher number state `\|20>`, and why the dimension `d` must exceed it |
27+| [`ex03_superposition_state.m`](ex03_superposition_state.m) | II.B | A normalized superposition of number states |
28+| [`ex04_coherent_state.m`](ex04_coherent_state.m) | II.C | A coherent state `\|alpha>` built by truncating its number-basis sum |
29+| [`ex05_thermal_state.m`](ex05_thermal_state.m) | II.D | A thermal state as a diagonal density matrix (mixed state) |
30+| [`ex06_squeezed_vacuum.m`](ex06_squeezed_vacuum.m) | II.E | A squeezed vacuum state — a superposition of even number states |
31+| [`ex07_nsfcs.m`](ex07_nsfcs.m) | II.F | A number-state-filtered coherent state (one Fock state removed) |
32+| [`ex08_atomic_states.m`](ex08_atomic_states.m) | II.G | Excited/ground states of a two-level atom and their superposition |
33+
34+## Operators (Sec. III)
35+
36+| Script | Article | What it shows |
37+| --- | --- | --- |
38+| [`ex09_field_operators.m`](ex09_field_operators.m) | III.A | Annihilation, creation and number operators in the number basis |
39+| [`ex10_field_hamiltonian.m`](ex10_field_hamiltonian.m) | III.B | The quantized electromagnetic-field Hamiltonian |
40+| [`ex11_pauli_operators.m`](ex11_pauli_operators.m) | III.C | Pauli `sigma_z`, raising/lowering operators, two-level-atom Hamiltonian |
41+
42+## Properties of quantum states (Sec. IV)
43+
44+| Script | Article | What it shows |
45+| --- | --- | --- |
46+| [`ex12_photon_distribution.m`](ex12_photon_distribution.m) | IV.A | Photon-number distributions `P_n` (bar plots): coherent vs. thermal |
47+| [`ex13_average_photons.m`](ex13_average_photons.m) | IV.B | Average photon number `<a'a>` for number, coherent, thermal states |
48+| [`ex14_g2_coherence.m`](ex14_g2_coherence.m) | IV.C | Second-order coherence `g^(2)(0)`: sub-Poissonian, Poissonian, super-Poissonian |
49+
50+## Atom–field interaction (Sec. V)
51+
52+| Script | Article | What it shows |
53+| --- | --- | --- |
54+| [`ex15_atom_field_rabi.m`](ex15_atom_field_rabi.m) | V | Jaynes–Cummings vacuum Rabi oscillations `P_e(t)`, `P_g(t)` |
55+| [`ex16_atomic_inversion.m`](ex16_atomic_inversion.m) | V | Atomic inversion `<sigma_z(t)>` with a coherent field: collapse and revival |
56+
57+## Two-mode field (Sec. VI)
58+
59+| Script | Article | What it shows |
60+| --- | --- | --- |
61+| [`ex17_coupled_cavities.m`](ex17_coupled_cavities.m) | VI.A | Photon exchange between two coupled cavities `P_{10}(t)`, `P_{01}(t)` |
62+
63+## A few numbl/MATLAB notes
64+
65+- **Time evolution** uses the matrix exponential `U = expm(-i*H*dt)`. The scripts
66+ build the propagator once, then apply it repeatedly in a loop.
67+- The examples set `hbar = 1` (and often `hbar*omega = 1`) so energies and times
68+ are in natural units, exactly as in the article.
69+- `i` and `j` both denote the imaginary unit. Scripts begin with `clear; clc;`
70+ following the article's convention.
71+- `ex16_atomic_inversion.m` is the heaviest example (a 100×100 propagator over a
72+ long evolution); give it a few seconds.
73+
74+## How it's deployed
75+
76+On every push to `main`, the workflow in
77+[`.github/workflows/deploy.yml`](.github/workflows/deploy.yml) bundles these
78+files with the numbl browser IDE and publishes the result to GitHub Pages — no
79+server, nothing to build by hand. Edit `numbl-project.json` to change the title
80+or which file opens first.
81+
82+## Reference
83+
84+Nilakantha Meher, *Quantum optics in MATLAB*, [arXiv:2309.14354](https://arxiv.org/abs/2309.14354).
ex01_number_states.madded+10−0View file
@@ -0,0 +1,10 @@
1+% Number states (Sec. II.A)
2+% A number state |n> is a column of the identity matrix: the (n+1)th column.
3+
4+clear; % Clear memory
5+clc; % Clear the command window/screen
6+d = 5; % dimension of the field
7+I = eye(d); % Identity matrix of dimension d
8+Vacuum = I(:,1) % First column of identity matrix: vacuum state |0>
9+Ket1 = I(:,2) % Second column of identity matrix: single-photon state |1>
10+Ket2 = I(:,3) % Third column of identity matrix: two-photon state |2>
ex02_number_state_20.madded+8−0View file
@@ -0,0 +1,8 @@
1+% Higher number state |20> (Sec. II.A)
2+% To represent |n>, the field dimension d must be larger than n.
3+
4+clear; % Clear memory
5+clc; % Clear the command window/screen
6+d = 21; % dimension of the field (must exceed 20)
7+I = eye(d); % Identity matrix of dimension d
8+Ket20 = I(:,21) % number state |20>: 21st column of the identity matrix
ex03_superposition_state.madded+14−0View file
@@ -0,0 +1,14 @@
1+% Superposition of number states (Sec. II.B)
2+% |psi> = (1/sqrt(3))|2> + (1/sqrt(2))|5> - (1/sqrt(6))|6>
3+% Any superposition of number states is a linear combination of identity
4+% columns; d must exceed the largest number state present.
5+
6+clear; % Clear memory
7+clc; % Clear the command window/screen
8+d = 7; % dimension of the field
9+I = eye(d); % Identity matrix of dimension d
10+Ket2 = I(:,3); % two-photon state |2>
11+Ket5 = I(:,6); % five-photon state |5>
12+Ket6 = I(:,7); % six-photon state |6>
13+Psi = 1/sqrt(3)*Ket2 + 1/sqrt(2)*Ket5 - 1/sqrt(6)*Ket6 % superposition state
14+N_psi = norm(Psi) % norm of the state Psi
ex04_coherent_state.madded+15−0View file
@@ -0,0 +1,15 @@
1+% Coherent state (Sec. II.C)
2+% |alpha> = exp(-|alpha|^2/2) * sum_n alpha^n/sqrt(n!) |n>
3+% The infinite sum is truncated at d-1; check that the norm stays ~1.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 10; % dimension of the field
8+I = eye(d);
9+alpha = 0.6; % amplitude of the coherent state
10+Coh = 0; % initialization
11+for x = 0:d-1
12+ Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I(:,x+1);
13+end
14+Coh % display the coherent state
15+N_c = norm(Coh) % check norm (should be very close to 1)
ex05_thermal_state.madded+15−0View file
@@ -0,0 +1,15 @@
1+% Thermal state (Sec. II.D)
2+% rho_th = 1/(1+nth) * sum_n (nth/(1+nth))^n |n><n|
3+% A mixed state: a square density matrix, diagonal in the number basis.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 10;
8+I = eye(d);
9+nth = 0.5; % average number of photons in the thermal state
10+RhoTh = 0; % initialization
11+for x = 0:d-1
12+ RhoTh = RhoTh + nth^(x)/(1+nth)^(x+1)*I(:,x+1)*I(:,x+1)';
13+end
14+RhoTh % display the thermal state
15+N_th = trace(RhoTh) % check trace (should be 1)
ex06_squeezed_vacuum.madded+17−0View file
@@ -0,0 +1,17 @@
1+% Squeezed vacuum state (Sec. II.E)
2+% |xi> = 1/sqrt(cosh r) * sum_n (-1)^n sqrt((2n)!)/(2^n n!) e^{i n theta} (tanh r)^n |2n>
3+% A superposition of the even number states only.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 20;
8+I = eye(d);
9+r = 0.3; % squeezing parameter
10+theta = pi/4; % squeezing direction
11+Sqz = 0; % initialization
12+for x = 0:(d/2)-1
13+ p = (1/sqrt(cosh(r)))*sqrt(prod(1:(2*x)))/(2^x*prod(1:x));
14+ Sqz = Sqz + p*(-1)^x*exp(i*x*theta)*(tanh(r))^x*I(:,2*x+1);
15+end
16+Sqz % squeezed vacuum state output
17+N_sqz = norm(Sqz) % check norm (should be 1)
ex07_nsfcs.madded+19−0View file
@@ -0,0 +1,19 @@
1+% Number state filtered coherent state, NSFCS (Sec. II.F)
2+% A coherent state with one number state |m> removed from the superposition,
3+% then renormalized.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 15;
8+m = 4; % number state |4> will be absent from the distribution
9+I = eye(d);
10+alpha = 0.8; % amplitude of the coherent state
11+nsfs = 0; % initialization
12+for x = 0:d-1
13+ if x == m
14+ nsfs = nsfs + 0*I(:,x+1);
15+ else
16+ nsfs = nsfs + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I(:,x+1);
17+ end
18+end
19+NSFS = nsfs/norm(nsfs) % normalize the state (note the 5th element is 0)
ex08_atomic_states.madded+11−0View file
@@ -0,0 +1,11 @@
1+% States of a two-level atom (Sec. II.G)
2+% Excited |e> and ground |g> states as 2-dimensional column vectors, and a
3+% normalized superposition |psi> = alpha|e> + beta|g>, |alpha|^2+|beta|^2 = 1.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+es = [1;0] % excited state |e>
8+gs = [0;1] % ground state |g>
9+alpha = sqrt(0.4); % superposition coefficient
10+beta = sqrt(0.6); % superposition coefficient
11+Psi = alpha*es + beta*gs % superposition state (norm = 1)
ex09_field_operators.madded+10−0View file
@@ -0,0 +1,10 @@
1+% Annihilation, creation and number operators (Sec. III.A)
2+% a|n> = sqrt(n)|n-1>, a'|n> = sqrt(n+1)|n+1>, N = a'a, N|n> = n|n>
3+% The creation operator is the adjoint of the annihilation operator.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 6; % dimension of the operator
8+A = diag(sqrt(1:d-1), 1) % Annihilation operator (superdiagonal)
9+Ad = A' % Creation operator (adjoint of A)
10+N = Ad*A % Number operator (diagonal: 0,1,2,...)
ex10_field_hamiltonian.madded+13−0View file
@@ -0,0 +1,13 @@
1+% Hamiltonian for the electromagnetic field (Sec. III.B)
2+% H = hbar*omega*(a'a + 1/2 I)
3+% We set hbar = omega = 1 to express energies in units of hbar*omega.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 6; % dimension of the operator
8+I = eye(d); % identity matrix
9+hbar = 1;
10+omega = 1; % frequency of the cavity field
11+A = diag(sqrt(1:d-1), 1); % Annihilation operator
12+Ad = A'; % Creation operator
13+H = hbar*omega*(Ad*A + (1/2)*I) % Hamiltonian (diagonal: 0.5, 1.5, 2.5, ...)
ex11_pauli_operators.madded+12−0View file
@@ -0,0 +1,12 @@
1+% Pauli operators and Hamiltonian for a two-level atom (Sec. III.C)
2+% H = (hbar*omega_0/2) sigma_z
3+% sigma_+ raises (|g> -> |e>) and sigma_- lowers (|e> -> |g>).
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+hbar = 1;
8+Omega_0 = 1; % atomic transition frequency
9+Sz = [1 0; 0 -1] % sigma_z operator
10+Splus = [0 1; 0 0] % raising operator sigma_+
11+Sminus = [0 0; 1 0] % lowering operator sigma_-
12+H = hbar*Omega_0/2*Sz % Hamiltonian for a two-level atom
ex12_photon_distribution.madded+38−0View file
@@ -0,0 +1,38 @@
1+% Photon number distribution (Sec. IV.A)
2+% P_n = |<n|psi>|^2. Poissonian for a coherent state, monotonic for a
3+% thermal state. Two bar plots are produced.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 15; % dimension of the field
8+I = eye(d);
9+
10+%%%%%%%% Coherent state
11+alpha = 2; % amplitude of the coherent state
12+Coh = 0; % initialization
13+for x = 0:d-1
14+ Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I(:,x+1);
15+end
16+for n = 0:d-1
17+ PnC(n+1) = norm(I(:,n+1)'*Coh)^2; % Probability of |n> in coherent state
18+end
19+
20+%%%%%%%% Thermal state
21+nth = 0.5; % average number of photons in the thermal state
22+RhoTh = 0;
23+for x = 0:d-1
24+ RhoTh = RhoTh + nth^(x)/(1+nth)^(x+1)*I(:,x+1)*I(:,x+1)';
25+end
26+for n = 0:d-1
27+ PnTh(n+1) = I(:,n+1)'*RhoTh*I(:,n+1); % Probability of |n> in thermal state
28+end
29+
30+%%%%%%%%%% Bar plots
31+n = 0:d-1;
32+figure(1)
33+bar(n, PnC)
34+title('Coherent state \alpha = 2'); xlabel('n'); ylabel('P_n')
35+
36+figure(2)
37+bar(n, PnTh)
38+title('Thermal state n_{th} = 0.5'); xlabel('n'); ylabel('P_n')
ex13_average_photons.madded+33−0View file
@@ -0,0 +1,33 @@
1+% Average number of photons (Sec. IV.B)
2+% <a'a> = <psi|a'a|psi> (pure states), Tr(a'a rho) (mixed states)
3+% Number state |4> -> 4, coherent (|alpha|^2 = 3) -> 3, thermal (nth) -> nth.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 20; % dimension of the field
8+I = eye(d); % identity matrix
9+
10+%%%% Number operator
11+A = diag(sqrt(1:d-1), 1); % Annihilation operator
12+Ad = A'; % Creation operator
13+AdA = Ad*A; % Number operator
14+
15+%% Number state |4>
16+Ket4 = I(:,5); % four-photon state |4>
17+AdANumber = Ket4'*AdA*Ket4 % average photon number in the number state
18+
19+%% Coherent state
20+alpha = sqrt(3); % amplitude of the coherent state
21+Coh = 0; % initialization
22+for x = 0:d-1
23+ Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I(:,x+1);
24+end
25+AdACoherent = Coh'*AdA*Coh % average photon number in coherent state (= |alpha|^2)
26+
27+%%% Thermal state
28+nth = 0.85; % assumed average number of photons in the thermal state
29+RhoTh = 0;
30+for x = 0:d-1
31+ RhoTh = RhoTh + nth^(x)/(1+nth)^(x+1)*I(:,x+1)*I(:,x+1)';
32+end
33+AdAThermal = trace(AdA*RhoTh) % average photon number in thermal state (= nth)
ex14_g2_coherence.madded+38−0View file
@@ -0,0 +1,38 @@
1+% Zero time-delay second-order coherence function g^(2)(0) (Sec. IV.C)
2+% g2(0) = <a'^2 a^2> / <a'a>^2
3+% Number state |n>: 1 - 1/n (sub-Poissonian); coherent: 1; thermal: 2.
4+
5+clear; % Clear memory
6+clc; % Clear the command window/screen
7+d = 25; % dimension of the field
8+I = eye(d); % identity matrix
9+A = diag(sqrt(1:d-1), 1); % Annihilation operator
10+Ad = A'; % Creation operator
11+AdA = Ad*A; % Number operator
12+Ad2A2 = Ad*Ad*A*A; % a'^2 a^2
13+
14+%% Number state |4>
15+Ket4 = I(:,5); % four-photon state |4>
16+AdANumber = Ket4'*AdA*Ket4;
17+Ad2A2Number = Ket4'*Ad2A2*Ket4;
18+g2Number = Ad2A2Number/(AdANumber)^2 % g2(0) for number state (= 1 - 1/4 = 0.75)
19+
20+%% Coherent state
21+alpha = sqrt(3); % amplitude of the coherent state
22+Coh = 0; % initialization
23+for x = 0:d-1
24+ Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I(:,x+1);
25+end
26+AdACoherent = Coh'*AdA*Coh;
27+Ad2A2Coherent = Coh'*Ad2A2*Coh;
28+g2Coherent = Ad2A2Coherent/(AdACoherent)^2 % g2(0) for coherent state (= 1)
29+
30+%%% Thermal state
31+nth = 0.85; % average number of photons in the thermal state
32+RhoTh = 0;
33+for x = 0:d-1
34+ RhoTh = RhoTh + nth^(x)/(1+nth)^(x+1)*I(:,x+1)*I(:,x+1)';
35+end
36+AdAThermal = trace(AdA*RhoTh);
37+Ad2A2Thermal = trace(Ad2A2*RhoTh);
38+g2Thermal = Ad2A2Thermal/(AdAThermal)^2 % g2(0) for thermal state (= 2)
ex15_atom_field_rabi.madded+41−0View file
@@ -0,0 +1,41 @@
1+% Atom-field interaction: vacuum Rabi oscillations (Sec. V)
2+% Jaynes-Cummings Hamiltonian; the atom starts excited with n photons in the
3+% cavity and periodically exchanges a photon with the field.
4+% Pe(t) = cos^2(g t sqrt(n+1)), Pg(t) = sin^2(g t sqrt(n+1))
5+
6+clear; % Clear memory
7+clc; % Clear the command window/screen
8+d = 10; % dimension of the cavity field
9+hbar = 1;
10+W0 = 1; % atomic frequency
11+Wf = 1; % cavity field frequency
12+g = 0.1; % coupling constant
13+A = diag(sqrt(1:d-1), 1); % Annihilation operator
14+Ad = A'; % Creation operator
15+Sz = [1,0;0,-1]; % sigma_z
16+Splus = [0,1;0,0]; % sigma_+
17+Sminus = [0,0;1,0]; % sigma_-
18+gs = [0;1]; % ground state
19+es = [1;0]; % excited state
20+I_a = eye(2); % identity operator for the atom
21+I_f = eye(d); % identity operator for the field
22+Hatom = (1/2)*hbar*W0*kron(Sz,I_f); % atomic Hamiltonian
23+Hfield = hbar*Wf*kron(I_a,Ad*A); % field Hamiltonian
24+Hint = hbar*g*(kron(Splus,A)+kron(Sminus,Ad)); % interaction Hamiltonian
25+H = Hatom + Hfield + Hint; % JC Hamiltonian
26+n = 4; % initial number of photons in the cavity
27+en = kron(es,I_f(:,n+1)); % atom excited, field has n photons |e,n>
28+gn = kron(gs,I_f(:,n+2)); % atom ground, field has n+1 photons |g,n+1>
29+Psi = en; % initial state
30+dt = 0.1; % time step
31+U = expm(-i*H*dt/hbar); % unitary time-evolution operator
32+T = 0:dt:30; % total evolution time
33+for t = 1:length(T)
34+ Pe(t) = norm(en'*Psi)^2; % probability of |e,n>
35+ Pg(t) = norm(gn'*Psi)^2; % probability of |g,n+1>
36+ Psi = U*Psi; % time-evolved state
37+ Psi = Psi/norm(Psi); % renormalize
38+end
39+plot(T, Pe, 'r', T, Pg, 'k') % Pe in red, Pg in black
40+title('Vacuum Rabi oscillations (n = 4, g = 0.1)')
41+xlabel('t'); ylabel('Probability'); legend('P_e', 'P_g')
ex16_atomic_inversion.madded+44−0View file
@@ -0,0 +1,44 @@
1+% Atomic inversion with a coherent field: collapse and revival (Sec. V)
2+% Atom excited, field in a coherent state. The atomic inversion <sigma_z(t)>
3+% collapses and later revives -- a hallmark of field quantization.
4+%
5+% Note: d = 50 and a long evolution make this the heaviest example here;
6+% give it a few seconds.
7+
8+clear; % Clear memory
9+clc; % Clear the command window/screen
10+d = 50; % dimension of the cavity field
11+hbar = 1;
12+W0 = 1; % atomic frequency
13+Wf = 1; % cavity field frequency
14+g = 0.1; % coupling constant
15+A = diag(sqrt(1:d-1), 1); % Annihilation operator
16+Ad = A'; % Creation operator
17+Sz = [1,0;0,-1]; % sigma_z
18+Splus = [0,1;0,0]; % sigma_+
19+Sminus = [0,0;1,0]; % sigma_-
20+gs = [0;1]; % ground state
21+es = [1;0]; % excited state
22+I_a = eye(2); % identity operator for the atom
23+I_f = eye(d); % identity operator for the field
24+Hatom = (1/2)*hbar*W0*kron(Sz,I_f);
25+Hfield = hbar*Wf*kron(I_a,Ad*A);
26+Hint = hbar*g*(kron(Splus,A)+kron(Sminus,Ad));
27+H = Hatom + Hfield + Hint; % JC Hamiltonian
28+alpha = 3; % coherent state amplitude
29+Coh = 0;
30+for x = 0:d-1
31+ Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I_f(:,x+1);
32+end
33+Psi = kron(es,Coh); % initial state: atom |e>, field coherent
34+dt = 0.1; % time step
35+U = expm(-j*H*dt); % unitary time-evolution operator
36+T = 0:dt:500; % total evolution time
37+for t = 1:length(T)
38+ W(t) = Psi'*kron(Sz,I_f)*Psi; % <sigma_z> (atomic inversion)
39+ Psi = U*Psi;
40+ Psi = Psi/norm(Psi);
41+end
42+plot(T, W)
43+title('Atomic inversion <\sigma_z(t)>: collapse and revival (\alpha = 3)')
44+xlabel('t'); ylabel('<\sigma_z>')
ex17_coupled_cavities.madded+31−0View file
@@ -0,0 +1,31 @@
1+% Coupled cavities: energy exchange between two field modes (Sec. VI.A)
2+% One photon starts in cavity 1, cavity 2 is empty. The photon is
3+% periodically transferred between the two cavities.
4+% P10(t) = cos^2(Jt), P01(t) = sin^2(Jt)
5+
6+clear; % Clear memory
7+clc; % Clear the command window/screen
8+d = 10; % dimension of the cavity field
9+hbar = 1;
10+W1 = 1; % resonance frequency of the first cavity
11+W2 = 1; % resonance frequency of the second cavity
12+J = 0.1; % inter-cavity coupling constant
13+A = diag(sqrt(1:d-1), 1); % Annihilation operator
14+Ad = A'; % Creation operator
15+I_f = eye(d); % identity operator for the fields
16+H = hbar*W1*kron(Ad*A,I_f) + hbar*W2*kron(I_f,Ad*A) + hbar*J*(kron(Ad,A)+kron(A,Ad));
17+S10 = kron(I_f(:,2),I_f(:,1)); % cavity 1 has |1>, cavity 2 has |0>
18+S01 = kron(I_f(:,1),I_f(:,2)); % cavity 1 has |0>, cavity 2 has |1>
19+Psi = S10; % initial state
20+dt = 0.1; % time step
21+U = expm(-i*H*dt/hbar); % unitary time-evolution operator
22+T = 0:dt:50; % total evolution time
23+for t = 1:length(T)
24+ P10(t) = norm(S10'*Psi)^2; % prob. photon in cavity 1
25+ P01(t) = norm(S01'*Psi)^2; % prob. photon in cavity 2
26+ Psi = U*Psi; % time evolution
27+ Psi = Psi/norm(Psi); % normalization
28+end
29+plot(T, P10, 'r', T, P01, 'k') % cavity 1 red, cavity 2 black
30+title('Photon exchange between coupled cavities (J = 0.1)')
31+xlabel('t'); ylabel('Probability'); legend('P_{10}', 'P_{01}')
numbl-project.jsonadded+4−0View file
@@ -0,0 +1,4 @@
1+{
2+ "title": "Quantum Optics in numbl",
3+ "entry": "README.md"
4+}