Quantum optics in MATLAB

We provide a MATLAB numerical guide at the beginner level to support students starting their research careers in theoretical quantum optics and related areas. These resources are also valuable for undergraduate and graduate students working on semester projects in similar fields.


Operations
nth column of a matrix A A(:,n) A⊗B (tensor product) kron(A,B) TABLE I. Used operations and their MATLAB commands.

I. INTRODUCTION
MATLAB is a user-friendly and robust framework for numerical computing based on matrix operations.Several numerical toolboxes or open-source packages written in MATLAB [1][2][3][4][5][6] have been designed to address analytically intractable problems in quantum mechanics, quantum optics and condensed matter physics.
In this tutorial, we present various numerical codes written in MATLAB to help students understand the basics of quantum optics.These codes can be easily extended to address a wide range of research problems in quantum optics and related areas that involve high-dimensional matrix manipulations.Importantlly, they can be executed in any version of MATLAB without requiring pre-installed packages.
Before we start writing codes in MATLAB, we follow the rules: • Pure quantum states are represented by normalized column matrices with unit norm, while mixed states are represented by square matrices with unit trace.
• Quantum operators are represented by square matrices.
• We need to pre-decide the dimension of the matrices for states and operators based on the specific problem at hand.The dimension refers to the number of rows or columns.If your results deviate from your expectations, you can consider increasing the dimension of the matrices.We will learn how to decide the appropriate dimensions of the matrices as we write code.
• We use 'clear' to clear memory and 'clc' to clear the output screen at the beginning of a code to prevent any numerical errors.
• To suppress output, add a semicolon ';' at the end of the respective line.
• The symbol % is used to comment a line.The sentence after % is used to explain the code to the readers.It will not be executed in the run of a program and will not be displayed at the output.
• We use "format short" for truncating the number of digits after the decimal point to be four.One may choose "formal long" to display more number of digits after the decimal point.
• We listed a few commands in Table .I that we use to perform calculations.

II. QUANTUM STATES IN MATLAB A. Number States
The number states are the most fundamental quantum states of a quantized electromagnetic (EM) field.The number state |n⟩ represents a field having n photons.
In MATLAB, we represent a number state |n⟩ with a column matrix in which the element at the (n + 1)th position is 1, while the other elements are 0. Accordingly, the vacuum state |0⟩ is represented by a column matrix in which the first element is 1 and the others are 0.
Trick: Consider each column of an identity matrix to represent a number state.The outputs are: It is to be noted that the above column matrices are derived from a 5×5 identity matrix (d=5).The output is a column matrix with 21st element is 1 and others are zero.The configuration represents the number state |20⟩.

B. Superposition of Number States
In the previous subsection, we mapped the number states to columns of an identity matrix.For example, n-photon state |n⟩ is mapped to (n + 1)th column of the identity matrix.Now, using this mapping, let us we write the following superposition state in MATLAB: We must ensure that the state |ψ⟩ is normalized.The norm of |ψ⟩ is calculated to be Note: Any superposition of number states can be expressed as a linear combination of the columns of an identity matrix.The dimension (d) of the matrix must be larger than the largest number-state present in the superposition state.

C. Coherent States
A coherent state |α⟩ is a superposition of all number states [7,8].In the number basis, it is expressed as [Ch. 3 of Ref. [9]] Therefore, to write this state in MATLAB, we need to sum all the columns of an identity matrix weighted by the coefficient e −|α| 2 /2 α n √ n! .In principle, one needs to add infinite number of such column matrices.However, for a given α, the coefficient e −|α| 2 /2 α n √ n! becomes negligible as n increases.Thus, we truncate the sum so the norm of the state remains very close to 1.Although we have truncated the above state up to d=10 (or n = 9), it is important to note that the norm N c of the coherent state is very close to 1.As we use "format short", the output displays N c=1.0000.If we use "format long", then it will display N c=0.999999999996367, which is very close to 1.

Code
If we take a larger value of 'alpha', the state with d=10 may not have the unit norm.For that case, we need to increase the value of 'd'.
Note: If the norm of a state is found to be less than 1, it is necessary to increase the dimension of the field 'd'.Please refer to A for an explanation of the significance of dimension 'd'.It illustrates the amount of error that may arise when choosing a smaller value for 'd'.

D. Thermal States
The electromagnetic radiation from an object at a non-zero temperature is a thermal light.The density matrix of a thermal state in the number basis is [Ch. 2 of Ref. [9]] where n th is the average number of photons in the thermal state ρ th .As thermal states are mixed, we represent them by a square matrices, and they are diagonal in number basis.As we can see from the above equation, it is a sum of square matrices |n⟩ ⟨n| for all n with appropriate coefficients (probabilities).To get the matrix form of |n⟩ ⟨n|, consider an example.For n = 3, we have It is a square matrix whose 4th element of the diagonal is 1.
Code for thermal state:

E. Squeezed Vacuum States
A squeezed vacuum state in the number basis is [Ch.7 of Ref. [9]] where r and θ are the squeeze parameters.This state is a superposition of all the even number states.

F. Number State Filtered Coherent State
A number state filtered coherent state (NSFCS) in number basis is [10] where the normalization constant

G. States of a two-level atom
Consider a two-level atom whose excited and ground states are represented by |e⟩ and |g⟩ respectively.We represent these states by column matrices of dimension 2. We have Then, any superposition of these states is The outputs are: We check that the norm of 'Psi' is 0.6325 2 + 0.7746 2 = 1.

III. OPERATORS IN MATLAB
The operators correspond to a quantized EM field will be written in number basis (|n⟩) and the operators correspond to two-level atoms will be written in atomic state basis (|e⟩ and |g⟩).

A. Annihilation, Creation and Number operators
The basic operators that we use to study the EM field in quantum optics are the annihilation operator, creation operator and the number operator [11].To represent these operators in matrix form, we first need to know their action on number states.
The action of annihilation and creation operators on a number state is The action of annihilation operator on a field state corresponds to the destruction of a photon from the field, while the action of creation operator corresponds to the creation of a photon in the field.The above two equations give the matrix form Note: The creation operator is adjoint of the annihilation operator.
The number operator is N = â † â, and its action on number state is Thus, it is diagonal in number basis.The expectation value of the number operator in a number state counts the number of photons in that state, that is, Code for annihilation, creation and number operators: The Hamiltonian for a quantized electromagnetic field is [Ch. 2 of Ref. [9], Ch. 1 of Ref. [11]] where ℏ = 6.626 × 10 −34 Js (Planck's constant) and ω is the frequency of the field.The operator Îf is an identity operator.In many textbooks, authors don't explicitly write the identity operator Îf .

Note:
We often take ℏω = 1 to normalize the output quantities in the unit of energy.
Code for EM field Hamiltonian: Let the excited and ground states of a two-level atom, represented by |e⟩ and |g⟩, have the energies ℏω 0 /2 and −ℏω 0 /2 respectively.Thus, the energy difference is ℏω 0 .The Hamiltonian for a two-level atom is [12] is a Pauli matrix.To account for the transition between the two energy levels, we consider the raising and lowering operators respectively.Their actions are Note: We often take ℏω 0 = 1 to normalize the output quantities in the unit of energy.
Code for atomic Hamiltonian, raising and lowering operators: Upon measurement on a number state |n⟩, we obtain n photons with unit probability.However, there are several nonclassical fields that do not possess a definite number of photons [13].When we measure the photon number in such fields, we collapse the field state to a particular number state with a definite number of photons.If we perform repeated measurements on an ensemble of such field states, we encounter different photon numbers with associated probabilities.For example, the photon-number probability distribution in a coherent state is Poissonian.
The probability of finding n photons in a given field state |ψ⟩ is In particular, the probability of detecting n photons in a coherent state |α⟩ is and in a thermal state is [Ch. 2 of Ref. [9]] Code for photon-number probability distributions in coherent and thermal states: The output bar plots are given in Fig. 1.The height of a bar indicates the probability of detecting a corresponding number state in coherent state (left) and thermal state (right).

B. Average number of photons
The average number of photons in a given field state |ψ⟩ is For number state |n⟩, which is the number of photons in the state.The average number of photons in a coherent state is [9] ⟨a and in a thermal state is Code for calculating the average number of photons: The outputs are: 3.0000 AdAThermal = 0.8500 In the above code, we selected the number state |4⟩ and calculated the average number of photons to be 4.This outcome illustrates that the number operator counts the number of photons in a number state.This is because the number states are eigenstates of the number operator.For coherent state, we take the amplitude to be α = √ 3 and calculated the average number of photons to be |α| 2 = 3.Similarly, we assumed n th to be 0.85 to create thermal state in the code and that is reflected at the output.

V. ATOM-FIELD INTERACTION
The strength of interaction (energy exchange) between an atom with an electromagnetic field in free space is small [15,16].This is enhanced by confining them in a cavity [Fig.2(a), [16,17]].The Hamiltonian for an atom-cavity system is [Ch. 4 of Ref. [9], Ch. 12 of Ref. [18]] This is the Jaynes-Cummings (JC) Hamiltonian [19].The first and second terms are the energy operators for a two-level atom and a cavity field, respectively.We have omitted the constant ℏω/2 from the energy operator of the cavity field since it doesn't significantly affect the time dynamics.The last term represents the interaction between field and atom, facilitating energy exchange in the form of photons.The atom absorbs a photon from the field and subsequently re-emits it into the field.
Important Note: When dealing with two physical systems in quantum mechanics, we need to be careful while writing their combined operators and states.An operator for a composite system is written as a tensor product of the operators of the individual sub-systems, that is, where Â and B are the operators corresponding to the sub-systems A and B, and Ĉ is an operator for the composite system.Similarly, a state of the combined system is also a tensor product of the states of the individual sub-systems, that is, It is important to note that we follow the order of the operators and the states to be the same, that is, we place the operator and state of the sub-system A before those of the sub-system B. This ensures that the operators of the subsystem A act meaningfully on the state of A, and similarly for the subsystem B. We may note the last term of Eq. ( 29), particularly the terms σ− â † and σ+ â, the atomic operators are placed first and then the field operators.To follow the same order, we re-write the JC Hamiltonian as Here, Îa and Îf are the identity operators correspond to two-level atom and the cavity field.The dimension of Îa is 2 and the dimension of Îf depends on the dimension of the field state.
To study the time dynamics of the system, we define the unitary operator to be where Ĥ is given in Eq. (29).
Let the combined initial state of atom and field be where the atom is in its excited state |e⟩ and the cavity field has precisely n photons.Note that we have followed the same order for the state too (atomic state and then field state), which is necessary to have a meaningful action of the operators on the state.For resonant case, that is, ω = ω 0 , the combined state evolves under the JC Hamiltonian as [Ch. 4 of Ref. [9]] where |e, n⟩ = |e⟩ ⊗ |n⟩ and |g, n + 1⟩ = |g⟩ ⊗ |n + 1⟩.At t = 0, the atom is in excited state and the cavity field is in n-photon state.At time t = π/(2g √ n + 1), the atom goes to its ground state by emitting a photon into the field.Then, at time t = π/(g √ n + 1), the atom comes back to its initial state |e⟩ by absorbing a photon from the field.Thus, the atom will be periodically exchanging energy with the field [see Fig. 2(b)].
The probability of finding the atom in excited state and the field in n-photon state is and the probability of finding the atom in ground state and the field in (n + 1)-photon state is Code for atom-field interaction:  The output plot is shown in Fig. 2(b).At t = 0, P e = 1 and P g = 0.
At t = π/(2g √ n + 1) = π/(2 × 0.1 × √ 4 + 1) ≈ 7, P g = 1 and P e = 0. Now, consider the initial state of the field to be in a coherent state and the atom in its excited state.We calculate the expectation value of the atomic operator σ z to be [Ch.4 of Ref. [9]] We consider the beam splitter to be 50:50 by taking θ = π/4.As a result, it splits the beam into two with equal intensities.
Let the input be |n, 0⟩, where n photons are in mode-a and the mode-b is in vacuum.Thus, the input intensity is ⟨â † â⟩ = n.The output state will be a superposition of all possible combinations generated from the distribution of n photons in two modes [Ch. 5 of Ref. [18]]: The average number of photons (∼ intensity) at the output port a out is and at the output port b out is After the 50:50 beam splitter transformation, each output gets half of its input intensity.
For n = 2, the input state is |2, 0⟩.The output state can be calculated using Eq. ( 46) to be Therefore, the probabilities of getting |2, 0⟩ , |1, 1⟩ and |0, 2⟩ are 1/4, 1/2 and 1/4, respectively.The average number of photons (∼ intensity) at the output port a out is and at the output port b out is For the input |n, m⟩, the output state after the beam splitter is [Ch. 5 of Ref. [18]] The outputs are: The Mach-Zehnder interferometer (MZI) is one of the simplest devices that demonstrates interference of an electromagnetic field at a single photon level [Ch.6 of [9]].It consists of two beam splitters, a few mirrors, and a phase shifter [See Fig. 4(b)].The first beam splitter divides the input field into two.A phase shifter in one of the arms/paths Quantum jump Unitary evolution FIG. 5.The probability of finding the atom in its excited state, calculated using master equation approach (black continuous line), as a function of time t in the presence of cavity dissipation.The result is compared with the Monte-Carlo wavefunction method (red dashed line).We set κ = 0.05, γ = 0, and g = 0.1.For MCWF, the number of realizations is N = 5000.The inset shows a single trajectory generated using MCWF method.A quantum jump occurs at t ∼ 56.time δt, the state |ψ(0)⟩ evolves under a non-Hermitian Hamiltonian where Ĥ is given in Eq. (32).Here, we consider only the cavity is dissipating.For δt ≪ 1, the state at a time δt is where Î = Îa ⊗ I f is the identity matrix for the atom-cavity system.As the Hamiltonian ĤNH is non-Hermitian, the evolution does not preserve the norm and the norm will be less than 1.So, the missing norm is which should be much less than 1.
• Once δp is calculated, it will be compared with a random number r chosen from a uniform distribution between 0 and 1.If r > δp, no jump occurs and the state |ψ(δt)⟩ will be normalized.And if r < δp, a quantum jump occurs, and the normalized state at δt will be • To get the state at 2δt, we need to follow the same procedure above: evolving under a non-Hermitian Hamiltonian, calculate the missing norm and compare with a random number to decide the state at a later time.This procedure has to be followed up to a time T = nδt.This forms a single trajectory.But, one has to get many such quantum trajectories from the initial state |ψ(0)⟩, and ensemble average over many such trajectories produces the exact evolution for the initial quantum state.
The following code considers the dissipation of an atom through cavity (κ ̸ = 0, γ = 0).The probability of finding the atom in excited state after each MCWF step is calculated through and taken averages over many realizations.The result is compared with the results of the master equation approach in Fig. 5.We take 5000 realizations.The inset shows one of the MCWF trajectories in which a quantum jump (photon loss) occurs at t ∼ 56.After the quantum jump, the atom and cavity reach to their ground state.If the results from MCWF method do not match with that of the master equation approach, one needs to increase the number of realizations and decrease the time step δt.However, for such cases, the MCWF method becomes computationally much more expensive than the Lindblad master equation method.

2 FIG. 2 .
FIG. 2. (a) An atom is trapped inside a cavity.(b) Time evolution of the probabilities Pe and Pg for atom-field coupling strength g = 0.1 and initial photon number n = 4. (c) Two cavities Cavity 1 and Cavity 2 are placed near to each other (coupled) to exchange energy.(d) Time evolution of the probabilities P10 and P01 for cavity-cavity coupling strength J = 0.1.The photon is periodically exchanged between the cavities.

For 50 :FIG. 4 .
FIG. 4. (a) A beam splitter: two input ports ain, bin and two output ports aout, bout.(b) A schematic of Mach-Zehnder interferometer consisting of two beam splitters, four mirrors and a phase-shifter.(c) Interference pattern: output intensity at aout as a function of ϕ.
Code for number states (vacuum state, single-photon state, and two-photon state): and hence, |ψ⟩ is normalized.As N psi=1, the state is normalized.The output state is a column matrix in which the 3rd, 6th and 7th elements are non-zero, which correspond to the amplitudes of the number states |2⟩, |5⟩, and |6⟩ in the superposition state, respectively.As the largest number-state present in the superposition state is |6⟩, the dimension (d) of the field is taken to be 7.In general, if a superposition state contains the largest number-state to be |n⟩, then 'd' must be larger than n+1.
2|α| 2m m! .This definition implies that the state is obtained if the number state |m⟩ is absent in the superposition.Note that the 5th element is zero indicating the absent of the number state |4⟩.
FIG.3.Time evolution of ⟨σz(t)⟩ (atomic inversion) when the atom is in excited state and field is in a coherent state.We set α = 3 and coupling strength g = 0.1.
Code for calculating expectation value of atomic operator σz when the field is in coherent state: