In [1]:

```
# All dependencies are included here
using PyPlot, NLsolve, LinearAlgebra, BenchmarkTools, SparseArrays, Arpack
# NLsolve is a key player in all of the routines here, and so I encourage you to read:
# https://github.com/JuliaNLSolvers/NLsolve.jl
```

In [2]:

```
# Fundamental constants (these are in "atomic units")
# https://en.wikipedia.org/wiki/Atomic_units
const c = 137.0; # speed of light
const ħ = 1.0; # planck constant
const ϵ0 = 1.0/4.0/pi; # permittivity of free space
const me = 1.0; # electron mass
const qe = 1.0; # electron charge
const nm = 18.904; # a nanometer in units of bohr radii
const μm = 18904.0; # a micrometer in units of bohr radii
const a0 = 5.29e-11; # bohr radius in SI units
const hartree_ev = 27.3114; # converts atomic unit energies to electron volts
const ev_hartree = 1.0/27.3114; # converts atomic unit energies to electron volts
const rads_ev = 1/(1.52e15); # converts frequency in radians/sec to electron volt energies
const ev_rads = 1.52e15; # converts electron volt energies to frequency in radians/sec
```

In [3]:

```
# The first struct constructs the 'electron class', working in the "velocity gauge" where
# momentum operators are important, while the second constructs the electron class in the
# "length gauge", where position operators are important.
struct electronic
N::Int64 # spectral dimension of the full hilbert space (vector space) describing the electron
Nred::Int64 # spectral dimension of truncated hilbert space (i.e., number of eigenvecs retained in computations)
H0::Matrix # matter Hamiltonian
p::Matrix # momentum operator
X::Array # a vector which holds all the eigenvectors and eigenvalues
end
struct electronic_r
N::Int64 # spectral dimension
Nred::Int64 # spectral dimension of truncated hilbert space (i.e., number of eigenvecs retained in computations)
H0::Matrix # matter Hamiltonian
r::Matrix # position
X::Array
end
# This struct below constructs the 'photon class', which refers to any kind of electromagnetic quanta essentially.
struct photonic
N::Int64 # spectral dimension
d::Int64 # polarization vector dimension
freqs::Vector # frequencies of modes 1 to N
modes::Array # photonic modes in absence of matter. rows are now components, columns are mode #
end
```

Given a spectral decomposition of the solutions of Maxwell's equations in terms of a complete and orthonormal set of vector modes $\mathbf{F}_n(\mathbf{r})$ with mode frequencies $\omega_n$, the spectral representation of the green's function is given by:

\begin{equation} \mathbf{G}(\mathbf{r},\mathbf{r}',\omega) = c^2\sum\limits_{n} \frac{\mathbf{F}^*(\mathbf{r})\otimes\mathbf{F}(\mathbf{r}')}{\omega_n^2-\omega^2-i\epsilon}, \end{equation}with $\epsilon$ a positive infinitesimal to enforce outgoing boundary conditions. These green's functions are used in the routine that calculates the response of the interacting electromagnetic field to a probe dipole.

In [12]:

```
function green!(i,j,ω,P::photonic,G=0.0)
# for now, positions pos1 and pos2 is going to be an integer index, which requires some input from the user
# for lower dimensional cases where G will be 1 or 2 components, just set other components to zero. That works.
spec_dim = P.N;
for counter = 1:spec_dim
G += conj(P.modes[i,counter])*(P.modes[j,counter])/((P.freqs[counter])^2-ω^2);
end
G = c^2*G;
return G
end
function dgreen!(i,j,ω,P::photonic,dG=0.0)
# for now, positions pos1 and pos2 is going to be an integer index, which requires some input from the user
# for lower dimensional cases where G will be 1 or 2 components, just set other components to zero. That works.
spec_dim = P.N;
for counter = 1:spec_dim
dG += -2*ω*conj(P.modes[i,counter])*(P.modes[j,counter])/((P.freqs[counter])^2-ω^2)^2;
end
dG = c^2*dG;
return dG
end
# a version that has explicit position dependence
function green2!(pos1,pos2,ω,P::photonic,G)
# for now, positions pos1 and pos2 is going to be an integer index, which requires some input from the user
# for lower dimensional cases where G will be 1 or 2 components, just set other components to zero. That works.
spec_dim = P.N;
space_dim = P.d;
for counter = 1:spec_dim
for ix = 1:space_dim
for iy = 1:space_dim
G[ix,iy] += conj(P.modes[pos1,counter,ix])*(P.modes[pos2,counter,iy])/((P.freqs[counter])^2-ω^2);
end
end
end
#G = c^2*G
return G
end
```

Out[12]:

This is used to solve the equation:

\begin{equation} \Big|I_3 + \lambda\mathbf{G}_n(\mathbf{r}_0,\mathbf{r}_0,\omega_n)\Big| = 0 \end{equation}where $\lambda$ is the coupling coefficient seen below.

In [5]:

```
function photon_potential_A2(str)
λ = str*qe^2/me/ϵ0/c^2;
return λ
end
```

Out[5]:

These functions calculate the "non-perturbative" interaction Hamiltonian (matrix) associated with virtual emission and absorption of photons by the electron. The interaction is calculated according to

\begin{equation} H_{\mathrm{int}}=\frac{q^2\hbar}{2m^2\epsilon_0}\sum\limits_{n=1}^{N_p}\sum\limits_{b=1}^{N_e}\frac{|\mathbf{F}_n(\mathbf{r}_0)|^2 p\psi_b\psi_b^{\dagger}p}{\omega_n(E_1-E_b-\omega_n)}, \end{equation}where the modes $\mathbf{F}_n$ (and their frequencies $\omega_n$) are stored in the photon class, while the vectors $\psi_b$, and the "energies" $E_b$ are eigenvectors and eigenvalues corresponding to:

\begin{equation} \left(H_0 + H_{\mathrm{int}}\right)\psi = E\psi. \end{equation}In practice, this nonlinear problem is solved as a fixed point problem by taking the eigenvectors and eigenvalues in $H_{\mathrm{int}}$ as specified by an initial guess. The momentum operator is given by $-i\hbar\nabla$, where $\nabla$ is the discretized gradient operator.

In [ ]:

```
# for now keeping things one dimensional in p dot A = px*Ax
function scf_potential2!(X,E::electronic,P::photonic,V_scf,str=1)
N = E.N;
H0 = E.H0;
p = E.p;
# X = E.X;
spec_dim = P.N;
vec_dim = P.d;
for count_el = 1:N
for count_pt = 1:spec_dim
tmp_factor = (ħ/2/ϵ0/(P.freqs[count_pt]))*abs(P.modes[1,count_pt])^2/(X[N^2+1]-X[N^2+count_el]-P.freqs[count_pt]);
V_scf += str*tmp_factor*p*(conj.(X[1+(count_el-1)*N:(count_el)*N]').*X[1+(count_el-1)*N:(count_el)*N])*p;
V_scf = real.(V_scf);
end
end
return V_scf
end
function scf_iter_eigs2(X,E::electronic,P::photonic,V_scf,str)
N = E.N;
H0 = E.H0;
p = E.p;
# X = E.X;
spec_dim = P.N;
vec_dim = P.d;
V_scf = scf_potential2!(X,E,P,V_scf,str);
V_scf = (V_scf+V_scf')/2;
H = H0 + V_scf;
vecs_H,vals_H = eigvecs(H),eigvals(H);
tmp1 = reshape(vecs_H,(N^2,1));
tmp2 = reshape(vals_H,(N,1));
X = [tmp1
tmp2];
# electron = electronic(N,H0,p,X)
return X
end
```

This program makes substantial use of routines for solving nonlinear equations: both nonlinear equations in a single variable, and nonlinear eigenvalue problems.In order to run the routines in this code, you must first Pkg.add("NLsolve").

Nonlinear equations in a single variable are solved with Newton's method. This is the default method.

Nonlinear eigenproblems are solved with fixed point iteration, which is achieved with fixedpoint(). The default method is method = :anderson with mixing parameter $m=0$ and damping parameter $\beta = 1$ (these parameters correspond to usual fixed-point iteration in which information about previous iterates is not explicitly used). A non-zero $m$ leads to a mixing of $m$ previous iterates in determining the next iterate.

Note that in some cases, we found that residuals of the first equation for the matter can vary wildly between iterations. In which case, consider using a non-zero $m$, as well as a damping parameter smaller than $1$.

In [15]:

```
function f_pt!(F,λ)#,str=1)
# weight = photon_potential_A2(str);
F[1] = -weight*green!(1,1,λ[1],pt_cav,0.0)[1] - 1
end
function j_pt!(J,λ)#,str=1)
# weight = photon_potential_A2(str);
J[1] = weight*dgreen!(1,1,λ[1],pt_cav,0.0)[1]
end
function construct_coupled_freqs(P::photonic,weight)
N = P.N;
freqs_int = zeros(N,1)
for count = 1:N
if abs(P.modes[count])^2 >= (1.0e-2)*abs(P.modes[1])^2
sol_pt = nlsolve(f_pt!,[(count+0.02)*pi*c/L_cav])
freqs_int[count] = sol_pt.zero[1];
else
freqs_int[count] = P.freqs[count]
end
end
freqs_int = reshape(freqs_int,(N));
freqs_int = sort(freqs_int);
return freqs_int
end
function mode(freq)
fieldmode = zeros(num_positions_cav);
for position = 1:num_positions_cav
fieldmode[position] = green2!(position,501,freq,pt_cav,[0.0])[1];
end
fieldmode = fieldmode / norm(fieldmode);
return fieldmode
end
function f_el2!(F,X)
X = scf_iter_eigs2(X,E,P,zeros(N_el,N_el),str);
# print("$(X)")
N = E.N;
# print("$(real(X[end-N+1:end])) \n")
for count = 1:size(X,1)
F[count] = X[count];#X[end-N+count]
end
end
```

Out[15]:

The routines below find the ground state of a strongly-coupled light-matter system in the non-perturbative regime by a fundamentally different method, called mean-field theory. In this scheme, the problem reduces like before to a nonlinear eigen-equation, but of a different form. Here, it is

\begin{equation} \left(-\frac{1}{2}\nabla^2 + U(\mathbf{r}) - q\mathbf{r}\cdot\langle \mathbf{E} \rangle \right)\psi = E\psi. \end{equation}where the expectation value of the quantized electric field is given by

\begin{equation} \langle \mathbf{E}(\mathbf{R}) \rangle = \sum\limits_n \sqrt{\frac{\hbar\omega_n}{2\epsilon_0}}\left(\alpha_n \mathbf{F}_n(\mathbf{R}) + c.c \right), \end{equation}where $\mathbf{R}$ is the center of charge of the system and the coefficients $\alpha_n$ are given by

\begin{equation} \alpha_n = q\frac{1}{2\epsilon_0\hbar\omega_n}\int d^3x ~ \mathbf{x}|\psi(\mathbf{x})|^2. \end{equation}As a result of the equations for $\psi$ featuring a term with $\psi(\mathbf{r})\int d^3x ~\mathbf{x}|\psi(\mathbf{x})|^2$, it is both nonlinear and non-local.

In [33]:

```
function scf_potential_meanfield!(X,E::electronic_r,P::photonic,V_scf,str=1)
# right now let's only worry about the ground state and see where that gets us
N = E.N;
Nred = E.Nred;
H0 = E.H0;
r = E.r;
# X = E.X;
spec_dim = P.N;
vec_dim = P.d;
Efield = 0;
d_av = qe*conj.(X[1:N]'*r*X[1:N]);
for count_pt = 1:spec_dim
Efield += 1/2/ϵ0*2*real(d_av)*abs(P.modes[1,count_pt])^2
end
V_scf = - qe*r*Efield
return V_scf
end
function scf_iter_eigs_meanfield(X,E::electronic_r,P::photonic,V_scf,str=1)
N = E.N;
Nred = E.Nred;
H0 = E.H0;
r = E.r;
# X = E.X;
spec_dim = P.N;
vec_dim = P.d;
V_scf = scf_potential_meanfield!(X,E,P,V_scf,str);
V_scf = (V_scf+V_scf')/2;
H = H0 + str*V_scf;
H = sparse(H)
if issparse(H) == true
vals_H, vecs_H = eigs(H,which=:SM,nev=Nred)
else
vals_H,vecs_H = eigvals(H),eigvecs(H);
end
#print("$(size(H))")
#print("$(size(vals_H))")
#print("$(size(vecs_H))")
tmp1 = reshape(vecs_H,(N*Nred,1));
tmp2 = reshape(vals_H,(Nred,1));
X = [tmp1
tmp2];
# electron = electronic_r(N,H0,r,X)
return X
end
function scf_photon_meanfield(E::electronic_r,P::photonic)
N = E.N;
Nred = E.Nred;
H0 = E.H0;
r = E.r;
X = E.X;
d_av = qe*(X[1:N]'*r*X[1:N]);
alpha = d_av/sqrt(2*ϵ0*ħ)*conj.(P.modes./sqrt.(P.freqs)');
# field_tmp = sqrt(ħ/2/ϵ0)*(sqrt.(P.freqs))(alpha.*(P.modes)+(conj.(alpha)).*conj.(P.modes))
# field = sum(field_tmp);
return alpha
end
function f_el_meanfield!(F,X)
N = E.N;
X = scf_iter_eigs_meanfield(X,E,P,zeros(N,N),str);
# print("$(X)")
Nred = E.Nred;
# print("$(real(X[end-Nred+1:end])) \n")
for count = 1:size(X,1)
F[count] = X[count];#X[end-N+count]
end
end
```

Out[33]:

Here, we show a few examples of electrons classes that can be generated. It is common in contexts of condensed matter physics to consider electrons specified by a discrete Hamiltonian in a "site" or "tight-binding" representation, in which the basis used to construct the Hamiltonian is a basis of states highly localized at a set of discrete lattice sites. The off-diagonal elements $H_{ij}$ determine the probability that an electron on site $j$ tunnels to site $i$. The on-diagonal elements $H_{ii}$ sometimes called on-site energies determine the energies of electrons at different sites in the absence of any tunneling. The functions below construct tight-binding electron models in both velocity and position gauges. More examples of electrons, not associated with any functions (just "manually built") are shown in the demo file.

In [9]:

```
function construct_matter(V,t,R,N)
Hm = SymTridiagonal(V,t*ones(N-1));
p = -1*im/R*SymTridiagonal(zeros(N),ones(N-1))
vals_Hm, vecs_Hm = eigvals(Hm), eigvecs(Hm);
tmp1 = reshape(vecs_Hm,(N^2,1));
tmp2 = reshape(vals_Hm,(N,1));
Xm = [tmp1
tmp2];
electron = electronic(N,N,Hm,p,Xm);
return electron
end
function construct_matter_r(V,t,R,N)
Hm = SymTridiagonal(V,t*ones(N-1));
positions = R*(-1:2/(N-1):1)
r = Diagonal(positions)
vals_Hm, vecs_Hm = eigvals(Hm), eigvecs(Hm);
tmp1 = reshape(vecs_Hm,(N^2,1));
tmp2 = reshape(vals_Hm,(N,1));
Xm = [tmp1
tmp2];
electron = electronic_r(N,N,Hm,r,Xm);
return electron
end
```

Out[9]:

Below we manually construct examples of photon systems which are quite relevant in current research:

- Photons in a metallic cavity which is long in an "axial" direction and small in two "transverse" directions.
- "Photons" (or electromagnetic fields) associated with collective oscillations of electrons in an electron gas (sometimes called plasmons).

For more physics details, see for example: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.80.245435

In [10]:

```
N_cav = 100;
d_cav = 1; # spatial dimension of the cavity mode, corresponds to number of vector components used to describe mode
L_cav = 18904.0; # 1 micron in atomic units
S_cav = (10.0*nm)^2; #area of cavity
emitter_pos = L_cav/2; # position where the emitter is located
freqs_cav = pi*(1:N_cav)*c/L_cav;
modes_cav = zeros(d_cav,N_cav);
for counter = 1:N_cav
modes_cav[1,counter] = sqrt(2.0/L_cav/S_cav)*sin(freqs_cav[counter]*emitter_pos/c); # for 2d or 3d, want two extra lines specifying those components
end
pt_cav = photonic(N_cav,d_cav,freqs_cav,modes_cav)
```

Out[10]:

In [11]:

```
# This would be an example of user input.
# All of these parameters are used to specify the spectral properties of the plasma waves
N_pl = 10000; # approximating a continuum of plasma waves
d_pl = 3; # spatial dimension of the plasma wave, corresponds to number of vector components used to describe mode
qmax_pl = 5.0/emitter_pos; # maximum spatial frequency of the plasma waves that I retain
β_pl = 1.0e10; # proportionality constant used to define the frequencies of the plasmons. SI units.
Nq_pl = 200; # sets number of spatial frequencies used to sample a continuum of spatial frequencies
Nθ_pl = 50; # sets number of propagation angles used to sample a continuum of propagation angles
qs_pl = 0:qmax_pl/Nq_pl:qmax_pl; #sets the spatial frequencies of the modes... what to do about angles??
θs_pl = 0:2*pi/Nθ_pl:2*pi; #sets their polar angle of propagation in the 2D plane
S_pl = (1.0/(qmax_pl/Nq_pl))^2;
ω_pl(q) = β_pl*sqrt(q/a0)*rads_ev*ev_hartree+1.0e-7;
index_pl(iq,iθ) = (Nθ_pl)*(iq-1) + iθ;
# All of these parameters are used to specify the spectral properties of the plasma waves
emitter_pos = 10*nm; # 10 nm away from the 2D plasma
# Creating all of the
freqs_pl = zeros(N_pl);
modes_pl = zeros(Complex{Float64},d_pl,N_pl);
for iq = 1:Nq_pl
q_tmp = qs_pl[iq];
norm_tmp = sqrt(q_tmp/2.0/S_pl)*exp(-q_tmp*emitter_pos);
for iθ = 1:Nθ_pl
θ_tmp = θs_pl[iθ];
counter = index_pl(iq,iθ)
freqs_pl[counter] = ω_pl(q_tmp);
modes_pl[:,counter] = norm_tmp*[cos(θ_tmp),sin(θ_tmp),im];
end
end
pt_pl = photonic(N_pl,d_pl,freqs_pl,modes_pl)
```

Out[11]:

- Rivera, Nicholas, Johannes Flick, and Prineha Narang. arXiv preprint arXiv:1810.09595 (2018). [Original impetus for this work.]

- Flick, Johannes, et al. Proceedings of the National Academy of Sciences 112.50 (2015): 15285-15290. [Another impetus for this work].

- Walker, Homer F., and Peng Ni. SIAM Journal on Numerical Analysis 49.4 (2011): 1715-1735. [Anderson acceleration and fixed point problems.]

- https://web.wpi.edu/Pubs/ETD/Available/etd-112409-140359/unrestricted/pni.pdf. [Fixed-point problems (and mixing) as used in quantum chemistry.]

- https://github.com/JuliaNLSolvers/NLsolve.jl. [NLsolve package and documentation.]