# Spin Torque Oscillations Triggered by In-plane Field

R. Arun1,R. Gopal2, V. K. Chandrasekar2, and M. Lakshmanan1 1Department of
Nonlinear Dynamics, School of Physics, Bharathidasan University,
Tiruchirapalli-620024, India
2Centre for Nonlinear Science & Engineering, School of Electrical &
Electronics Engineering, SASTRA Deemed University, Thanjavur- 613 401, India.

###### Abstract

We study the dynamics of a spin torque nano oscillator that consists of
parallelly magnetized free and pinned layers by numerically solving the
associated Landau-Lifshitz-Gilbert-Slonczewski equation in the presence of a
field-like torque. We observe that an in-plane magnetic field which is applied
for a short interval of time ($<$1ns) triggers the magnetization to exhibit
self-oscillations from low energy initial magnetization state. Also, we
confirm that the frequency of oscillations can be tuned over the range
$\sim$25 GHz to $\sim$72 GHz by current, even in the absence of field-like
torque. We find the frequency enhancement up to 10 GHz by the presence of
field-like torque. We determine the Q-factor for different frequencies and
show that it increases with frequency. Our analysis with thermal noise
confirms that the system is stable against thermal noise and the dynamics is
not altered appreciably by it.

##  I Introduction

The spin torque nano oscillator (STNO) is a well-known candidate for its
potential application for microwave generation in the GHz range Grollier ;
Slavin ; Locatelli . Basically, the nanosized STNO consists of two
ferromagnetic layers, a free layer, and pinned layer, where the former one can
change its direction of magnetization while the latter one corresponds to a
fixed magnetization. The free layer is comparatively thinner than the pinned
layer. These two layers are separated by a nonmagnetic but conductive or
insulative layer called a spacer. The magnetization of the free layer can be
made to oscillate by tuning the current and/or magnetic field and the
magnetization oscillations can be converted into voltage oscillations by the
Giant Magnetoresistive effect Slon ; Berger ; Kiselev ; Rippard .

In an STNO the magnetization of the free layer is manipulated by the transfer
of spin angular momentum between the incoming spin polarized electrons and the
local magnetic moments in the free layer Slon ; Berger . The spin-torque
induced precession of the free layer’s magnetization has been experimentally
verified in the presence of current and magnetic field up to 26 GHz Kiselev ;
Mancoff1 ; Houss ; Mancoff2  and around 18 GHz in magnetic tunnel junctions
Nazarov1 ; Nazarov2 . The theoretical investigation on STNO predicts that the
spin torque acts like an external power against damping of the magnetization
and maintains self-oscillations. In the case of self-oscillations, the
magnetization never converges to a final state as if the system had no damping
at all. With appropriate conditions between the spin torque, magnetic field
and damping constant, it is possible to form a very stable precession even in
the presence of large thermal fluctuations Li .

The STNOs are characterized by various parameters such as frequency, power,
and Q-factor. Their efficient performance arises essentially due to their
nano-scale dimension and ability to tune their oscillation frequencies with
the enhancement of power and Q-factor. Apart from the requirement of high
frequency, tunability for a large range of frequencies by current or field is
also desired for different applications Z_Zeng .

Investigations on STNOs with perpendicularly magnetized free layer and
parallelly magnetized pinned layer in the presence of an out-of-plane field
show the possibility of (1) steady state precession with a maximum frequency
of 6.3 GHz with a high power of 0.55 $\mu$W Kubota  and Q-factor 135 and (2)
steady state precession with a maximum frequency of 45 GHz with a high
Q-factor over 3000 and power 64 nW Maehara . Without an out-of-plane field,
the steady state precession is not possible for this configuration of STNO,
and it can be induced by appropriate field-like torque Taniguchi1 ; Guo .

A further enhancement of frequency in the STNOs with parallelly magnetized
free layer has been achieved by applying an in-plane field Bonetti ; Suzuki1 ;
Suzuki2 ; ZM_Zeng ; Arun1  or by tilting the pinned layer’s magnetization to
out-of-plane Zhou1 ; Zhou2 ; Zhou3 ; Zhou4 ; Lv ; Arun2 . In particular,
Bonetti $et~{}al$. have experimentally observed the tunability of the
frequency by current, magnetic field strength and magnetic field direction in
nanocontact based STNOs Bonetti . They have found frequencies up to 46 GHz for
the strengths of magnetic field and current, ranging from 8 kOe to 14 kOe and
from 10 mA to 20 mA, respectively. Arun $et~{}al$., have studied the STNO in
the presence of an in-plane magnetic field and found high frequency
oscillations up to 68 GHz Arun1 . Also, these authors have identified the
oscillations with frequency up to 75 GHz in the absence of the field by
tilting the polarization of the pinned layer to out-of-plane Arun2 .

Therefore, the high frequency self-oscillations up to or above 70 GHz in an
STNO are achieved only by applying continuous in-plane/out-of-plane field
Bonetti ; Arun1  or by tilting the polarization of the pinned layer Arun2 .
Especially, getting self-oscillations corresponding to the low energy initial
magnetization state (i.e. aligned along easy axis direction) is practically a
challenging task. On the contrary, achieving the self-oscillations by current
without applying continuous external field and tilting the polarization of the
pinned layer still remains to be investigated for the case of parallelly
magnetized free layer. Therefore, in this present study, we examine the
occurrence of self-oscillations in an STNO having parallelly magnetized free
and pinned layers in the presence of a field-like torque. We show the rapid
response in the magnetization of the free layer towards the transition from a
steady state (no oscillatory state) to a self-oscillatory state after being
triggered by the in-plane magnetic field applied for a short interval of time.
The steady state is an equilibrium state where the magnetization converges to
settle finally. We find that the in-plane field applied for a short duration
causes the onset of self-oscillations of the magnetization with frequency up
to 72 GHz, even in the absence of field-like torque. Also, we show the
stability of the STNO against thermal noise and increment of Q-factor with
frequency.

The paper is organized as follows. Section II addresses the geometry of the
STNO and the governing equation, namely the, Landau-Lifshitz-Gilbert-
Slonczewski (LLGS) equation along with the details of the effective field. The
occurrence of self-oscillations of the magnetization due to the in-plane
field, the tunability of the frequency by current and field-like torque are
presented in Section III. Finally, the conclusions are summarized in Section
IV. We present the salient features of the derivation of the LLGS equation in
Appendix A and then in Appendix B we provide the impact of the thermal noise
on the oscillation frequency and power spectral density.

##  II Model Description and Governing Equation

Figure 1: Schematic picture of the model

We consider an STNO as shown in Fig. 1 with free and pinned layers (mild gray
color), and a spacer (dark gray color). The unit vectors of the magnetization
of the pinned layer and free layer are denoted by ${\bf p}$ and ${\bf m}$,
respectively.

The LLGS equation that governs the dynamics of the magnetization of the free
layer is given by (See Appendix A),

$\displaystyle\frac{d{\bf m}}{dt}=-\gamma~{}{\bf m}\times{\bf
H}_{eff}+\alpha~{}{\bf m}\times\frac{d{\bf m}}{dt}+\gamma H_{S}~{}{\bf
m}\times({\bf m}\times{\bf p})+\gamma\beta H_{S}~{}{\bf m}\times{\bf p}.$ (1)

In Eq. (1), ${\bf m}=m_{x}~{}{\bf e}_{x}+m_{y}~{}{\bf e}_{y}+m_{z}~{}{\bf
e}_{z}$ or ${\bf m}=(m_{x},m_{y},m_{z})$, $|{\bf m}|$ = 1, where ${\bf
e}_{x}$, ${\bf e}_{y}$ and ${\bf e}_{z}$ are the unit vectors along the
positive $x$, $y$ and $z$ directions, respectively. $\gamma$ is the
gyromagnetic ratio, $\alpha$ is the Gilbert damping parameter, $H_{S}$ is the
strength of the spin-transfer torque and $\beta$ is the ratio between the
strengths of field-like torque Li ; Zhang  and the spin-transfer torque. The
sign and magnitude of the field-like torque depend on the material parameters
such as diffusion constant, spin-flip relaxation time, thickness of the free
layer and the exchange interaction between the itinerant electron and the
magnetic background Zhang ; Shpiro ; Li1 . $H_{S}$ is defined as

$\displaystyle H_{S}=\frac{H_{S0}}{1+\lambda~{}{\bf m.p}},$ (2)

where $H_{S0}=\hbar\eta I/2eM_{s}V$. Here $I$ is the current passing through
the free layer, $\hbar(=h/2\pi)$ is the reduced Plank’s constant, $e$ is the
electron charge, $M_{s}$ is the saturation magnetization of the free and
pinned layers, $V$ is the volume of the free layer, $\eta$ and $\lambda$ are
dimensionless parameters that determine the magnitude and the angular
dependence of the spin-transfer torque, respectively. Here, the pinned layer’s
magnetization is fixed along the negative $x$-direction, and therefore ${\bf
p}=-{\bf e}_{x}$. The positive (negative) current corresponds to the flow of
electrons from the free (pinned) to pinned (free) layer. The effective field
${\bf H}_{eff}$ of the free layer is given by

$\displaystyle{\bf H}_{eff}~{}=~{}$ $\displaystyle H_{k}m_{x}~{}{\bf
e}_{x}-4\pi M_{s}m_{z}~{}{\bf e}_{z}+H_{a}(t)~{}{\bf n},$ (3)

where $H_{k}$ is the magneto-crystalline anisotropy field, $4\pi M_{s}$ is the
demagnetization field and $H_{a}(t)$ is the externally applied time-dependent
in-plane field, which can be applied continuously or for a short interval of
time. The latter case may be defined as a short in-plane field since it is
applied for a few nanoseconds only and it triggers self-oscillations as we see
in the next section. The time up to which the short in-plane field is applied
can be called cut-off time and it is denoted by $T_{off}$. In Eq. (3), ${\bf
n}=\cos\phi_{H}~{}{\bf e}_{x}+\sin\phi_{H}~{}{\bf e}_{y}$ is the unit vector
along which the in-plane field $H_{a}$ is applied and $\phi_{H}$ is the in-
plane field angle between ${\bf n}$ and ${\bf e}_{x}$ Taniguchi2 . The
effective field is defined by ${\bf H}_{eff}=-\partial E/\partial(M_{s}{\bf
m})$, where the $E$ is the effective energy density given by Taniguchi3

$\displaystyle E=$ $\displaystyle-\frac{M_{s}}{2}[H_{k}({\bf m}.{\bf
e}_{x})^{2}-4\pi M_{s}({\bf m}.{\bf e}_{z})^{2}]-M_{s}H_{a}{\bf m}.{\bf n}.$
(4)

In our present study, we consider the situation where the pinned and free
layers are made up of cobalt and the material parameters corresponding to
cobalt are chosen as Taniguchi2  $M_{s}$ = 1448.3 emu/c.c., $H_{k}$ = 18.6
kOe, $\eta$ = 0.54, $\lambda=\eta^{2}$, $\gamma$ = 17.64 Mrad/(Oe s), $\alpha$
= 0.005, $\mu_{0}$ = 1 and $V$ = $\pi\times 60\times 60\times 2$ nm3.

##  III Results

###  III.1 Oscillations triggered by in-plane magnetic field

Figure 2:  Temporal evolution of $m_{x}$ (a) in the absence of in-plane field,
(b) in the presence of in-plane field. Temporal evolution of (c) $m_{x}$ and
(d) energy density in the presence of short in-plane field with the cut-off
time $T_{off}$ = 1 ns, $H_{a}$ = 10 kOe and the direction $\phi_{H}$ = 30∘.
The mild and dark gray background colors in (c) indicate the cut-off time
$T_{off}$ up to which the in-plane field is applied and the minimum cut-off
time $T^{min}_{off}$, respectively. Here $\beta$ = 0.5 and $I$ = -15.5 mA.

Equation (4) implies that there exist two low energy magnetization states
${\bf m}_{1}$ = (1, 0, 0) and ${\bf m}_{2}$ = (-1, 0, 0) (we call as primary
equilibrium states) in the absence of current and in-plane field, marked by
the black solid points in Figs.3(a) and 3(b). Here, ${\bf m}_{2}$ is the
linearly stable equilibrium state, and ${\bf m}_{1}$ is the unstable
equilibrium point in the presence of current alone (i.e. in the presence of
current and absence of the external field). It means that if we choose the
initial magnetization state nearer to these states the system asymptotically
reaches the equilibrium state ${\bf m}_{2}$ (Fig. 2(a)). Though the state
${\bf m}_{2}$ is linearly stable it is nonlinearly unstable, which means that
if we choose the initial magnetization state far away from ${\bf m}_{1}$ or
${\bf m}_{2}$ , the system settles into a new stable self-oscillatory state in
the asymptotic limit. However, for $I\neq$0, $H_{a}\neq$0 and $\phi_{H}\neq$0,
both ${\bf m}_{1}$ and ${\bf m}_{2}$ are no longer equilibrium states for the
system. In this case, the system exhibits a new equilibrium state ${\bf
m}_{3}$ (and we call this as the secondary equilibrium state), which is stable
and the value of the equilibrium point depends on the value of the current and
applied field (See Fig. 2(b)). Interestingly, we observe that the system shows
oscillations for the initial magnetization state chosen nearer to ${\bf
m}_{1}$ (or ${\bf m}_{2}$) by applying the external magnetic field for a short
interval of time.

Figure 3: (a) Magnetization trajectories corresponding to the magnetization
reversal from $m_{x}$ = +1 to -1 and (b) basins of self-oscillation state
(gray points) and steady state (white region) in the absence of field when $I$
= -15.5 mA and $\beta$ = 0.5. Magnetization trajectories of (c) the steady
state due to the continuously applied in-plane field for $H_{a}$ = 10 kOe and
$\phi_{H}$ = 30∘ and (d) a trajectory of the self-oscillation due to the short
in-plane field with $H_{a}$ = 10 kOe, $\phi_{H}$ = 30∘ and $T_{off}$ = 1 ns.
Here $\beta$ = 0.5. The dotted black circles in (b) and (c) correspond to
$\theta=\theta_{osc}$.

In order to demonstrate the above, we fix the values of current the and the
field-like torque as $I=-15.5$ mA, $\beta$=0.5 respectively, and examine the
dynamics of the system in the absence and presence of the magnetic field in
Figs.2 and 3. For $H_{a}=0$, when the initial magnetization is chosen near the
unstable fixed point ${\bf m}_{1}$, the system shows magnetization reversal
and asymptotically reaches the stable state ${\bf m}_{2}$ as shown in
Figs.2(a) and 3(a). Fig. 3(b) shows the basin stability of the fixed point
${\bf m}_{2}$ as a white region, that is if the initial magnetization state is
taken in this region, the magnetization reaches the state ${\bf m}_{2}$
asymptotically. If we choose the initial magnetization state in the gray
region which refers to the basin of self-oscillatory state in Fig. 3(b), the
system will show self-oscillations of the magnetization. When we apply the in-
plane field continuously with $H_{a}$ = 10 kOe and, $\phi_{H}=30^{o}$, where
the initial magnetization is again chosen near to the unstable fixed point
${\bf m}_{1}$, the system reaches a secondary equilibrium state ${\bf
m}_{3}=(-0.22,-0.23,0.95)$ asymptotically, which is clearly evident from Fig.
3(c). The corresponding time evolution of $m_{x}$ and $m_{z}$ are plotted in
Figs.3(c) and 4, respectively.

However, if we apply a short duration in-plane field, for example in the form

$\displaystyle{\bf H_{a}}(t)=\begin{cases}10(\cos 30^{o}~{}{\bf e}_{x}+\sin
30^{o}~{}{\bf e}_{y})\text{ kOe},&\text{ $t\leq T_{off}$}\\\ 0,&\text{
$t>T_{off}$}\end{cases}$

the system shows self-oscillations even when the initial magnetization is
chosen in the basin of ${\bf m}_{2}$. The magnetization trajectory
corresponding to the self-oscillations of the magnetization is plotted after
leaving out the transients in Fig. 3(d) for $T_{off}=1~{}ns$.

The mechanism behind the emergence of self-oscillations of the magnetization
is as follows. When we apply the in-plane field the system reaches the
secondary equilibrium state and it emerges in the basin of self-oscillatory
state. So, when we cut off the field the system shows self-oscillations even
for the initial magnetization state chosen in the basin of ${\bf m}_{2}$.

Now, the temporal evolution of $m_{z}$ corresponding to the steady state
approaching towards the equilibrium state ${\bf m}_{3}$, which is also plotted
in Fig. 4, confirms that the magnetization settles into ${\bf m}_{3}$ after 2
ns. This implies that the system will exhibit self-oscillations on applying
the short duration in-plane field with the cut-off time $T_{off}\geq$ 2 ns.
Since the self-oscillations is achieved only after $T_{off}$, reducing the
cut-off time $T_{off}$ will reduce the delay in getting oscillations.
Therefore it is essential to find out the minimum cut-off time
$T^{min}_{off}$, so that the self-oscillations can be achieved with
$T_{off}\geq T^{min}_{off}$.

Figure 4: The temporal evolution of $m_{z}$ in the presence of continuous in-
plane field. The horizontal dashed line corresponds to
$m_{z}=cos(\theta_{osc})$ for finding the minimum cut-off time
$T_{off}^{min}$. The gray background indicates the time up to $T_{off}^{min}$.

The $T^{min}_{off}$ is the time at which the magnetization enters the basin of
self-oscillatory state and continues to precess inside the basin, around the
new steady state ${\bf m}_{3}$. One of the ways to estimate $T^{min}_{off}$ is
by identifying the zenith angle $\theta_{osc}$ (see Fig. 3(b)) in such a way
that the region corresponding to 0${}^{\circ}\leq\theta\leq\theta_{osc}$ does
not include the basin of steady state. For instance, the dotted black circle
corresponding to $\theta_{osc}$ = 39∘ drawn in Fig. 3(b) does not include the
basin of steady state. Here, $T^{min}_{off}$ can be considered as the time at
which the magnetization enters the dotted circle and continues to precess
inside the circle, around ${\bf m}_{3}$. The time at which the magnetization
enters the dotted circle is the same as the time at which the $m_{z}$ begins
to oscillate above the horizontal line drawn corresponding to
$m_{z}=\cos\theta_{osc}$ (see Fig. 4). For instance, the minimum cut-off time
$T^{min}_{off}$ is identified as 0.57 ns corresponding to the $\theta_{osc}$ =
39∘ as shown in Fig. 4.

Figure 5: The minimum cut-off time $T_{off}^{min}$ against (a) $\beta$ for
$H_{a}$ = 10 kOe and (b) $H_{a}$ for $\beta$ = 0.5. Here, $I$ = -15.5 mA.

The oscillations of the magnetization for $T_{off}$ = 1 ns are confirmed in
Fig. 2(c), where $m_{x}$ is plotted for $I$ = -15.5 mA, $\beta$ = 0.5,
$H_{a}(t)$ = 10 kOe($t\leq T_{off}$) and 0 Oe($t>T_{off}$), and $\phi_{H}$ =
30∘. We can observe that the oscillations are damped out as $t\rightarrow
T_{off}$ and triggered back as $t>T_{off}$ due to the short in-plane field. By
using the same procedure followed to determine the $T_{off}^{min}$ for $\beta$
= 0.5 we can also identify the values of $T_{off}^{min}$ for different values
of $\beta$ and $H_{a}$ and these are plotted in Fig.5(a) and Fig.5(b),
respectively. Figs.5 show that the minimum cut-off time does not exceed
$\sim$0.6 ns, which implies that the in-plane field need not be applied beyond
0.6 ns to trigger the self-oscillations. The slight change in $T_{off}^{min}$
due to $\beta$ and $H_{a}$ can be attributed to the fact that the size of the
basin of self-oscillatory state does not change considerably due to the change
in $\beta$ or $H_{a}$.

The time evolution of the energy density during the process of oscillations
due to the short in-plane field shown in Fig. 2(c) is plotted in Fig. 2(d).
Fig. 2(d) implies that in the presence of field and current the energy
oscillates initially and then tends to reach steady state energy of 0.8 J/cm3.
Consequently, the magnetization damps to reach the steady state ${\bf m}_{3}$
within the basin of the self-oscillatory state as shown in Figs.3(c) and 4.
After the in-plane field is switched off at time $T_{off}$(=1 ns) the energy
density emerges back to steadily oscillate between 1.50 J/cm3 and -0.88 J/cm3
and therefore $m_{x}$ oscillates steadily after 1 ns as shown in Fig.2(c).

Figure 6: Oscillation frequency with respect to (a) $\beta$ and $\phi_{H}$
when $T_{off}$ = 1 ns and $I$ = -15.5 mA and (b) $\beta$ and $I$ when $H_{a}$
= 10 kOe, $\phi_{H}$ = 30∘ and $T_{off}$ = 1 ns. The numerically (solid lines)
and analytically (open circles) obtained oscillation frequency in the presence
of short in-plane field against (c) current and (d) $\beta$ for $H_{a}$ = 10
kOe, $\phi_{H}$ = 30∘ and $T_{off}$ = 1 ns. Power spectral density in the
presence of short in-plane field for different values of (e) $\beta$ while $I$
= -15.5 mA and (f) current while $\beta$ = 0 for $\phi_{H}$ = 30∘. Here,
$T_{off}$ = 1 ns.

###  III.2 Frequency and power spectral density

We also find that the self-oscillations of the magnetization are not triggered
for all the directions of the short duration in-plane field. This is because
the position of the steady state ${\bf m}_{3}$ formed by the in-plane field
depends on the direction of the short duration in-plane field. For some
directions of the short duration in-plane field, the state ${\bf m}_{3}$ is
formed inside the basin of the steady state and, consequently, no oscillation
is exhibited. In order to verify the range of $\phi_{H}$ for which the self-
oscillation is possible, we plot the oscillation frequency between $\beta$ and
$\phi_{H}$ in Fig. 6(a) for $T_{off}$ = 1 ns and $I$ = -15.5 mA. The values of
($\beta,\phi_{H}$) corresponding to the white region between $\phi_{H}$ = 90∘
and 270∘ exhibit self-oscillations and the magnetization settles into ${\bf
m}_{2}$. This is because when $\phi_{H}$ is between 90∘ and 270∘, ${\bf
m}_{3}$ is formed outside the basin of self-oscillatory state and therefore no
self-oscillatory state is identified in the presence of the short in-plane
field. The system exhibits oscillations for $\phi_{H}$ between 0∘ and 90∘ and
270∘ and 360∘, where it is observed that the frequency is independent of the
in-plane field angle $\phi_{H}$.

Similarly, the current dependence of the self-oscillations is also identified
by plotting the frequency against $\beta$ and $I$ for $\phi_{H}$ = 30∘ and
$T_{off}$ = 1 ns in Fig. 6(b) and it is observed that the oscillatory region
lies between the nonoscillatory regions (white regions). From Fig. 6(b) we can
see that the frequency is minimum at $I=I_{c}^{low}$ and maximum at
$I=I_{c}^{upp}$. For $\beta=0$, $I_{c}^{low}$ is -12.4 mA and $I_{c}^{upp}$ is
-17.3 mA. Thus we conclude that self-oscillations are possible only when the
magnitude of the current is above the lower critical current ($I_{c}^{low}$)
and below the upper critical current ($I_{c}^{upp}$). The values of both the
critical currents decrease when the field-like torque is negative ($\beta<0$).
Further, from Fig. 6(b) we can also observe an enhancement in the frequency as
a function of the current and field-like torque. The frequency spectra for
different values of field-like torque and current are plotted in Fig. 6(c) and
(d) against the current and $\beta$, respectively. In Fig. 6(c) we have
plotted the numerically computed frequency against the current by solid lines
for different values of $\beta$ = 0, 0.5 and -0.5 in the presence of the short
duration in-plane field with $\phi_{H}$ = 30∘ and $T_{off}$ = 1 ns. From this
figure, we can confirm a large range of tunability in frequency from $\sim$25
GHz to $\sim$72 GHz by the current, and this is important for applications for
producing a wide range of frequencies. The open circles correspond to the
frequencies obtained from the expression

$\displaystyle f=\frac{\gamma\widetilde{m}_{z}}{2\pi(1+\alpha^{2})}$
$\displaystyle\left[-\frac{H_{k}}{2}-4\pi
M_{s}+\frac{H_{S0}(\alpha-\beta)}{\lambda(1-\widetilde{m}_{z}^{2})}\left(\frac{1}{\sqrt{1-\lambda^{2}(1-\widetilde{m}_{z}^{2})}}-1\right)\right],$
(5)

solved with the approximation $dm_{z}/dt=0$ by the procedure followed in
Ref.[31]. In the expression (5), $\widetilde{m}_{z}$ is the average value of
$m_{z}$ computed over $n$ number of precessions after reaching the self-
oscillations. The small deviation between the analytical and numerical
frequencies may be attributed to the approximation $m_{z}(t)=constant$ used in
this analysis. Fig. 6(c) implies that the self-oscillations of the
magnetization can be triggered by the short duration in-plane field even in
the absence of field-like torque and that $\beta<0$ is advantageous over
$\beta>0$ for getting the high frequency oscillations. To confirm the
occurrence of self-oscillations for the appropriate range of $\beta$ (from
-0.5 to 0.5, for Cobalt), we have plotted the graph of frequency against
$\beta$ for different values of current, namely $I$ = -15 mA, -15.5 mA and -16
mA, in Fig. 6(d) for the short in-plane field with $\phi_{H}$ = 30∘ and
$T_{off}$ = 1 ns. As we can see in Fig. 6(d), for a given current, the field-
like torque can increase the frequency up to 10 GHz.

The Q-factor (ratio between the peak frequency and line-width) of the STNO is
expected to be large for the betterment of its applications. In order to
determine the Q-factor for different frequencies, we have plotted the power
spectral density corresponding to different values of field-like torque for
$I$ = -15.5 mA and current for $\beta=0$ in Figs.6(e) and (f), respectively,
in the presence of the short in-plane field with $\phi_{H}$ = 30∘ and
$T_{off}$ = 1 ns. From the power spectral density plotted in Figs.6(e), the
Q-factor for the peaks corresponding to $\beta$ = 0, 0.5 and -0.5 are
estimated to be 471.62, 636.82 and 753.11, respectively. Similarly, from Fig.
6(f), the Q-factor for the peaks corresponding to $I$ = -13.5 mA, -15.0 mA and
-16.5 mA are given by 463.92, 591.63 and 743.42, respectively. The estimated
values of the Q-factor imply that it gets enhanced considerably with the
enhancement of the frequency. Also in Appendix B the impact of the thermal
noise on the oscillation frequencies and the power spectral density is
discussed.

The role of the in-plane magneto-crystalline anisotropy in the free layer is
important for the oscillations triggered by the short in-plane field. We have
verified that when the magneto-crystalline anisotropy is along the positive
z-direction the short-inplane field triggers no magnetization oscillations and
the magnetization only settles at (1,0,0) for positive current and at (-1,0,0)
for negative current irrespective of field-like torque. The free layer with
out-of-plane magneto-crystalline anisotropy exhibits the possibility of
oscillations with a frequency range of $\sim$6 GHz only in the presence of a
continuous perpendicular field Arun1 .

###  III.3 Impact of free layer volume on critical currents

The theory discussed above focuses on the spin-torque oscillator dynamics of
the magnetization vector ${\bf m}$, corresponding to the volume of the free
layer $V$ = 2$\times\pi\times$60$\times$60 nm3. We have found that the
magnitude of the critical currents $I_{c}^{low}$ and $I_{c}^{upp}$ are 12.4 mA
and 17.3 mA, respectively, in the presence of the short duration in-plane
field with $\phi_{H}$ = 30∘ and $T_{off}$ = 1 ns for $\beta$ = 0. The above
range of current (above 10 mA) and field (10 kOe) may be considered to be
quite high but we wish to note here that such values have been operated for
STNO previously. Bonetti $et~{}al$., have experimentally observed such
oscillations in STNO with the in-plane field 14.5 kOe and current 20 mA
Bonetti . Also, Taniguchi $et~{}al$., have theoretically studied for the
relaxation time in STNO for the range of current 0 to 25 mA Taniguchi2 .
Nevertheless, the magnitudes of the critical currents can be decreased further
by decreasing the volume of the free layer. To confirm that this is indeed
true, we have plotted the critical currents $I_{c}^{low}$ and $I_{c}^{upp}$
for different volumes of the free layer, namely $V_{1}$ =
2$\times\pi\times$60$\times$60 nm3 Arun1 ; Arun2 ; Taniguchi2 , $V_{2}$ =
2$\times\pi\times$50$\times$50 nm3, $V_{3}$ = 1$\times\pi\times$60$\times$60
nm3, $V_{4}$ = 2.5$\times$64$\times$64 nm3 Li  and $V_{5}$ =
3$\times$30$\times$30 nm3 Sbiaa . The results are shown in Fig. 7 which
essentially imply that the critical currents at which the system exhibits
self-oscillations of the magnetization that are directly proportional to the
volume of the free layer. Thus, a decrease in the volume of the free layer
decreases the critical currents required for the onset and offset of the self-
oscillations.

Figure 7: Critical currents for the onset and offset of the steady state
precession for different volumes of the free layer in the presence of the
short duration in-plane field. Here $H_{a}$ = 10 kOe, $\phi_{H}$ = 30∘,
$T_{off}$ = 1 ns and $\beta$ = 0.

##  IV Conclusions

We have studied an STNO with parallelly magnetized free and pinned layers in
the presence of field-like torque by solving the LLGS equation. We have
systematically shown the methodology to obtain self-oscillations of the free
layer’s magnetization by applying an in-plane field for a short duration of
time ($<$1ns) corresponding to the low energy initial magnetization state.
Expression for the frequency was derived and its validity was confirmed by
numerical simulation. We may emphasize that the self-oscillations of the
magnetization manifest even in the absence of field-like torque and one can
realize the possibility of enhancement in the oscillation frequency up to 10
GHz in the presence of field-like torque. We have also explored the tunability
of the frequency from $\sim$25 GHz to $\sim$72 GHz by the current. We
determined the Q-factor for different frequencies and confirmed that the
Q-factor increases with frequency. Finally, we have also confirmed that
thermal noise does not affect the dynamics and the oscillation frequency
appreciably.

## Acknowledgements

The works of V.K.C. and R. G are supported by the DST-SERB-CRG Grant No.
CRG/2020/004353 and they wish to thank DST, New Delhi for computational
facilities under the DST-FIST programme (SR/FST/PS-1/2020/135) to the
Department of Physics. M.L. wishes to thank the Department of Science and
Technology for the award of DST-SERB National Science Chair.

## Appendix A: Landau-Lifshitz-Gilbert-Slonczewski equation

In this Appendix, we present some brief details on the derivation of the LLGS
equation. For more details, we may refer to Lakshman . The magnetic moment of
an electron due to its orbital motion ${\bm{\mathcal{M}}_{orb}}$ in a circular
orbit is given by

$\displaystyle{\bm{\mathcal{M}}_{orb}}=I_{e}{\bf A},$ (A.1)

where $I_{e}$ is the current due to the orbiting electron and ${\bf A}$ is the
vector area of the orbital in which the electron is revolving. For an
electron, ${\bf A}=-\pi r^{2}\hat{L}$, where $r$ is radius of the circular
orbit and $\hat{L}$ is the unit vector of orbital angular momentum of the
electron. If the electron moves with a uniform speed $v$ in the circular orbit
with a time period $T$, we can write $I_{e}$ as, $I_{e}=e/T={e}/{(2\pi r/v)}$.
By substituting $I_{e}$ and ${\bf A}$ in Eq.(A.1), we can get

$\displaystyle{\bm{\mathcal{M}}_{orb}}=-\frac{evr}{2}~{}\hat{L}.$ (A.2)

For the circular orbit motion, the orbital angular momentum of the electron L
is given by ${\bf L}=rmv~{}\hat{L}$, where $m$ is mass of the electron. Now,
Eq.(A.2) can be written as ${\bm{\mathcal{M}}_{orb}}=-({\mu_{B}}/{\hbar}){\bf
L}$, where ${\mu}_{B}=e\hbar/2m$ is the Bohr magnetron. Analogously, one can
define the magnetic moment due to the spin of an electron as

$\displaystyle{\bm{\mathcal{M}}_{spin}}=-\gamma^{\prime}~{}{\bf S},$ (A.3)

where $\gamma^{\prime}=g\frac{\mu_{B}}{\hbar}$. $\gamma^{\prime}$ is the
gyromagnetic ratio, $g$ is the g-factor and ${\bf S}$ is the spin angular
momentum. The torque ${\bf N}$ exerted on the spin magnetic moment which is
immersed in the magnetic field $\bm{{H}}$ is given by

$\displaystyle{\bf N}=\frac{d{\bf
S}}{dt}=\mu_{0}~{}{\bm{\mathcal{M}}_{spin}}\times{{\bf H}},$ (A.4)

where $\mu_{0}$ is permeability in vacuum. By substituting ${\bf S}$ from
Eq.(A.3) in Eq.(A.4) we can derive

$\displaystyle\frac{d{\bm{\mathcal{M}}_{spin}}}{dt}=-\gamma^{\prime\prime}~{}{\bm{\mathcal{M}}_{spin}}\times{{\bf
H}},$ (A.5)

where $\gamma^{\prime\prime}=\gamma^{\prime}\mu_{0}$. For each spin magnetic
moment within a small volume $dV$ the Eq.(A.5) can be written as

$\displaystyle\frac{d{\bm{\mathcal{M}}_{spin,i}}}{dt}=-\gamma^{\prime\prime}~{}{\bm{\mathcal{M}}_{spin,i}}\times{{\bf
H}}.$ (A.6)

Here the small volume $dV$ is chosen in such a way that the magnetic field is
uniform inside it. If we take volume average on both sides of Eq.(A.6)

$\displaystyle\frac{d{\left(\sum_{i}\bm{\mathcal{M}}_{spin,i}/dV\right)}}{dt}=-\gamma^{\prime\prime}~{}{\left(\sum_{i}\bm{\mathcal{M}}_{spin,i}/dV\right)}\times{\bf{H}}.$
(A.7)

The quantity $\sum_{i}\bm{\mathcal{M}}_{spin,i}/dV$ is the net magnetic moment
per unit volume and it is denoted as the magnetization ${\bf M}$. Then
Eq.(A.7) can be rewritten as

$\displaystyle\frac{d\bf M}{dt}=-\gamma^{\prime\prime}~{}{\bf M}\times{\bf
H}.$ (A.8)

Here ${\bf M}=M_{s}{\bf m}$, where $M_{s}$ is the magnitude of the
magnetization called saturation magnetization and ${\bf m}$ is the unit
magnetization vector. Equation (A.8) represents the steady precession of ${\bf
M}$ around the field ${\bm{H}}$. But in reality the magnetization will damp
and align with the field after finite time due to energy dissipation. This
damping can be represented phenomenologically by introducing an additional
term in Eq.(A.8) as Ralph ; Landau ; Gilbert

$\displaystyle\frac{d\bf M}{dt}=-\gamma^{\prime\prime}~{}{\bf M}\times{\bf
H}-\frac{\lambda}{M_{s}}~{}{\bf M}\times\left({\bf M}\times{\bf H}\right).$
(A.9)

Equation (A.9) is called Landau-Lifshitz equation was first derived by Landau
and Lifshitz in 1935 to understand the dynamics of the magnetization in the
presence of magnetic field for small damping Landau ; Lakshman . The damping
term in Eq.(A.9) was replaced with a more convincing form for large damping in
1955 by Gilbert Gilbert  and the modified Landau-Lifshitz equation is given by

$\displaystyle\frac{d\bf M}{dt}=-\gamma~{}{\bf M}\times{\bf
H}+\frac{\alpha}{M_{s}}~{}{\bf M}\times\frac{d{\bf M}}{dt}.$ (A.10)

The two forms Eq.(A.9) and Eq.(A.10) are known to be equivalent with
$\gamma^{\prime\prime}=\gamma/(1+\alpha^{2})$ and
$\lambda=\gamma\alpha/(1+\alpha^{2})$. After dividing by $M_{s}$, Eq. (A.10)
is rewritten for the evolution of ${\bf m}$ as

$\displaystyle\frac{d\bf m}{dt}=-\gamma~{}{\bf m}\times{\bf H}+\alpha~{}{\bf
m}\times\frac{d{\bf m}}{dt}.$ (A.11)

Eq. (A.10) or Eq. (A.11) is now called the Landau-Lifshitz-Gilbert equation
and the magnetic field ${\bf H}$ is now replaced by the ${\bf H}_{eff}$ which
includes the magneto-crystalline anisotropy field, demagnetization field and
external field for a uniformly magnetized material. Later, in order to study
the influence of current polarized along with the direction ${\bf p}$ on the
magnetization an in-plane spin-transfer torque term $\gamma H_{S}~{}{\bf
m}\times({\bf m}\times{\bf p})$ was introduced by Slonczewski in 1996 Slon .
Then to investigate the impact of a perpendicular spin-transfer torque the
term $\gamma\beta H_{S}~{}{\bf m}\times{\bf p}$ was introduced by Zhang
$et~{}al.,$ in 2002 Zhang ; Li1 . With all these torques, the Landau-Lifshitz-
Gilbert-Slonczewski equation is given by

$\displaystyle\frac{d{\bf m}}{dt}=-\gamma~{}{\bf m}\times{\bf
H}_{eff}+\alpha~{}{\bf m}\times\frac{d{\bf m}}{dt}+\gamma H_{S}~{}{\bf
m}\times({\bf m}\times{\bf p})+\gamma\beta H_{S}~{}{\bf m}\times{\bf p}.$
(A.12)

This is the starting equation (1) for the unit vector of the magnetization.

## Appendix B: Impact of Thermal Noise

Figure 8: (a) Oscillation frequency due to short duration in-plane field
corresponding to the temperatures $T$ = 0 K (solid line) and $T$ = 300 K (open
circle). (b) Power spectral density in the presence of thermal noise
corresponding to $T$ = 300 K and $I$ = -15.5 mA. Here $T_{off}$ = 1 ns,
$H_{a}$ = 10 kOe and $\phi_{H}$ = 30∘.

For practical applications, the dynamics of an STNO are required to be stable
against the temperature. Therefore, it is important to verify the impact of
the thermal noise on the frequency of the magnetization oscillations and this
is studied by including the thermal field due to the thermal noise in the
effective field as follows:

$\displaystyle{\bf H}_{eff}~{}=~{}$ $\displaystyle H_{k}m_{x}~{}{\bf
e}_{x}-4\pi M_{s}m_{z}~{}{\bf e}_{z}$
$\displaystyle+H_{a}(t)(\cos\phi_{H}~{}{\bf e}_{x}+\sin\phi_{H}~{}{\bf
e}_{y})+{\bf H}_{th},$ (B.1)

where the thermal field is given by

$\displaystyle{\bf H}_{th}=\sqrt{D}~{}{\bf G},~{}~{}~{}~{}~{}D=\frac{2\alpha
k_{B}T}{\gamma M_{s}\mu_{0}V\triangle t}.$ (B.2)

In the above, ${\bf G}$ is the Gaussian random number generator vector of the
oscillator with components $(G_{x},G_{y},G_{z})$, which satisfies the
statistical properties $<G_{m}(t)>=0$ and
$<G_{m}(t)G_{n}(t^{\prime})>=\delta_{mn}\delta(t-t^{\prime})$ for all
$m,n=x,y,z$. $k_{B}$ is the Boltzmann constant, $T$ is the temperature,
$\mu_{0}$ is the magnetic permeability in free space and $\triangle t$ is the
step size of the time scale used in the simulation.

To figure out the impact of the thermal noise on the frequency and power
spectral density, Figs.8(a) and (b) are plotted for different values of field-
like torques in the presence of short duration in-plane field corresponding to
$T_{off}$ = 1 ns and $\phi_{H}$ = 30∘ for the temperatures $T$ = 0 K and $T$ =
300 K. In Fig.8(a), the frequencies corresponding to $T$ = 0 K and 300 K have
been represented by the solid lines and the open circles, respectively. From
Fig.8(a) and by comparing Figs.8(b) and 6(e), it can be understood that the
results are not altered by the thermal noise in any perceptible way.

## References

  * (1) J. Grollier, V. Cros and A. Fert, Phys. Rev. B 73, 060409(R) (2006).
  * (2) A. Slavin and V. Tiberkevich, IEEE. T. Magn. 45, 1875 (2009).
  * (3) N. Locatelli and V. Cros, J. Grollier, Nat. Mater. 13, 11 (2014).
  * (4) J.C. Slonczewski, Phys. Rev. B 39 6995 (1989); J. Magn. Magn.Mater. 159, (1996).
  * (5) L. Berger, Phys. Rev. B 54, 9353 (1996).
  * (6) S. I. Kiselev, J. C. Sankey, I. N. Krivorotov, N. C. Emley, R. J. Schoelkopf, R. A. Buhrman, and D. C. Ralph, Nature London 425, 380 (2003).
  * (7) W. Rippard, M. Pufall, and A. Kos, Appl. Phys. Lett. 103, 182403 (2013).
  * (8) F. B. Mancoff, N. D. Rizzo, B. N. Engel, and S. Tehrani, Nature 437, 393 (2005).
  * (9) D. Houssameddine, U. Ebels, B. Delaet, B. Rodmacq, I. Firastrau, F. Ponthenier, M. Brunet, C. Thirion, J.-P.Michel, L. Prejbeanu-Buda, M.-C. Cyrille, O. Redon, and B. Dieny, Nat. Mater. 6, 447 (2007).
  * (10) F. B. Mancoff, N. D. Rizzo, B. N. Engel, and S. Tehrani, Appl. Phys. Lett. 88, 112507 (2006).
  * (11) A. V. Nazarov, H. M. Olson, H. Cho, K. Nikolaev, Z. Gao, S. Stokes, and B. B. Pant, Appl. Phys. Lett. 88, 162504 (2006).
  * (12) A. V. Nazarov, K. Nikolaev, Z. Gao, H. Cho, and D. Song, J. Appl. Phys. 103, 07A503 (2008).
  * (13) Z. Li and S. Zhang, Phys. Rev. B 68, 024404 (2003).
  * (14) Z. Zeng, G. Finocchio, B. Zhang P. Khalili Amiri, J. A. Katine, I. N. Krivorotov, Y. Huai, J. Langer, B. Azzerboni, K. L. Wang, and H. Jiang, Sci. Rep. 3, 1426 (2013).
  * (15) H. Kubota, K. Yakushiji, A. Fukushima, S. Tamaru, M. Konoto, T. Nozaki, S. Ishibashi, T. Saruya, S. Yuasa, T. Taniguchi, H. Arai, and H. Imamura, Appl. Phys. Express 6, 103003 (2013).
  * (16) H. Maehara, H. Kubota, Y. Suzuki, T. Seki, K. Nishimura, Y. Nagamine, K. Tsunekawa, A. Fukushima, H. Arai, T. Taniguchi, H. Imamura, K. Ando, and S. Yuasa, Appl. Phys. Express 7, 023003 (2014).
  * (17) T. Taniguchi, S. Tsunegi, H. Kubota, and H. Imamura, Appl. Phys. Lett. 104, 152411 (2014).
  * (18) Y.-Y. Guo, H.-B. Xue, and Z.-J. Liu, AIP Advances 5, 057114 (2015).
  * (19) S. Bonetti, P. Mudali, F. Manoff, J. Akerman, Appl. Phys. Lett. 94, 102507 (2009).
  * (20) H. Suzuki, H. Endo, T. Nakamura, T. Tanaka , M. Doi, S. Hashimoto, H. N. Fuke, M. Takagishi, H. Iwasaki, M. Sahashi, J. Appl. Phys. 105, 07D124 (2009).
  * (21) H. Suzuki, T. Nakamura, H. Endo, M. Doi, H. Tsukahara, H. Imamura, H. N. Fukem, S. Hashimoto, H. Iwasaki, M. Sahashi, Appl. Phys. Lett. 99, 092507 (2011).
  * (22) Z. M. Zeng, P. Upadhyaya, P. Khalili Amiri, K. H. Cheung, J. A. Katine, J. Langer, K. L. Wang, and H. W. Jiang, Appl. Phys. Lett. 99, 32503 (2011).
  * (23) R. Arun, R. Gopal, V. K. Chandrasekar, and M. Lakshmanan, J. Magn. Magn. Mater. 532, 167989 (2021).
  * (24) Y. Zhou, C. L. Zha, S. Bonetti, J. Persson, and J. Akerman, Appl. Phys. Lett. 92, 262508 (2008).
  * (25) Y. Zhou, S. Bonetti, C. L. Zha, and J. Akerman, New J. Phys. 11, 103028 (2009).
  * (26) Y. Zhou, C. L. Zha, S. Bonetti, J. Persson, and J. Akerman, J. Appl. Phys. 105, 07D116 (2009).
  * (27) Y. Zhou, H. Zhang, Y. Liu, and J. Aherman, J. Appl. Phys. 112, 063903 (2012).
  * (28) G. Lv, H. Zhang, X. Cao, Y. Liu, Z. Hou, Y. Qin, G. Li, and L. Wang, AIP Adv. 5, 077171 (2015).
  * (29) R. Arun, R. Gopal, V. K. Chandrasekar, and M. Lakshmanan, J. Appl. Phys. 127, 153903 (2020).
  * (30) M. Lakshmanan, Phil. Trans. R. Soc. A 369, 1280 (2011).
  * (31) D. C. Ralph and M. D. Stiles, J. Magn. Magn. Mater. 320, 1190 (2008).
  * (32) L. D. Landau and E. M. Lifshitz, Phys. Z. Sowjet., bf 8, 153 (1935).
  * (33) T. L. Gilbert, IEEE Trans. Magn. 40, 3443 (2004).
  * (34) S. Zhang, P. M. Levy and A. Fert, Phys. Rev. Lett. 88, 236601 (2002).
  * (35) A. Shpiro, P. M. Levy, and S. Zhang, Phys. Rev. B 67, 104430 (2003).
  * (36) Z. Li, S. Zhang, Z. Diao, Y. Ding, X. Tang, D. M. Apalkov, Z. Yang, K. Kawabata, and Y. Huai, Phys. Rev. Lett. 100, 246602 (2008).
  * (37) T. Taniguchi, H. Kubota, Phys. Rev. B 93, 174401 (2016).
  * (38) T. Taniguchi, S. Tsunegi, H. Kubota, and H. Imamura, J. Appl. Phys. 117, 17C504 (2015).
  * (39) T. Taniguchi, T. Ito, S. Tsunegi, H. Kubota, and Y. Utsumi, Phys. Rev. B 96, 024406 (2017).
  * (40) R. Sbiaa, Curr. Appl. Phys. 14, 1521 (2014).


# Time resolved optical response of the Dicke’s model via the nonequilibrium
Green’s function approach

###### Abstract

Due to their conceptual appeal and computational convenience, two-level
systems (TLS) and their generalisations are often used to investigate
nonlinear behavior in quantum optics, and to assess the applicability of
theoretical methods. Here the focus is on second harmonic generation (SHG)
and, as system of interest, on the Dicke model, which consists of several TLSs
inside an optical cavity. The main aspect addressed is the scope of non-
equilibrium Green’s function (NEGF) to describe the effect of disorder and
electron-electron (e-e) interactions on the SHG signal. For benchmarking
purposes, exact diagonalization (ED) results are also presented and discussed.
SHG spectra obtained with NEGF and ED are found to be in very good mutual
agreement in most situations. Furthermore, inhomogeneity in the TLS and e-e
interactions reduce the strength of SHG, and the reduction is stronger with
inhomogeneity than with interactions. This trend is consistently noted across
different (small to large) system sizes. Finally, a modified NEGF approach is
proposed to account for cavity leakage, where the quantum photon fields are
coupled to a bath of classical oscillators. As to be expected, within this
mixed quantum-classical scheme a decrease in the intensity of the fluorescent
spectra takes place depending on the entity of cavity leakage.

###### keywords:

second harmonic generation, Dicke model, Disordered Dicke model, Non-
equilibrium Green’s functions, Exact diagonalization, Cavity leakage

Megha Gopalakrishna*, Yaroslav Pavlyukh*, Claudio Verdozzi*

Megha Gopalakrishna
Department of Physics, Division of Mathematical Physics,
Lund University, 22100 Lund, Sweden
Email Address<EMAIL_ADDRESS>

Yaroslav Pavlyukh
Institute of Theoretical Physics, Faculty of Fundamental Problems of
Technology,
Wroclaw University of Science and Technology, 50-370 Wroclaw, Poland
Email Address<EMAIL_ADDRESS>

Claudio Verdozzi
Department of Physics, Division of Mathematical Physics and ETSF,
Lund University, 22100 Lund, Sweden
Email Address<EMAIL_ADDRESS>

##  1 Introduction

Second harmonic generation (SHG) is a non-linear phenomenon in which two
incident photons are absorbed simultaneously, and a single photon is emitted
with doubled frequency [1]. SHG is the focus of extensive basic and applied
research in physics, chemistry, biology, engineering and medicine [2, 3, 4, 5,
6], because of its fundamental interest as optical process, and because its
vast use in characterization technique and optical devices. Extensive
investigations of SHG have been conducted both theoretically and
experimentally in many systems and there is now a rather comprehensive and
firm understanding of this optical process. There are, however, still some
particular situations where our grasp of SHG remains conceptually not
complete. The one of interest in this work is when many photon effects play a
role, and yet quantum fluctuations remain noticeable [7, 8, 9]. This paper
focusses on this specific regime, using the Dicke model [10] as the system of
interest.

The Dicke model holds significance in cavity optics, and has been vastly used
in the study of superradiance [10, 11]. It is also a go-to model for many
theoretical explorations concerning basic aspects of light-matter interaction,
and has also been considered in a simplified form by applying the so-called
rotating wave approximation (RWA)[12, 13]: in this case, the model is mostly
known as as the Tavis-Cummings model [14] (see e.g. [15] for a review).

In more detail, the original Dicke model [10, 11] describes a cavity mode
interacting with a series of identical two-level systems (TLSs), and the
original Dicke Hamiltonian is

$\displaystyle\hat{H}_{DM}=\omega_{a}\hat{a}^{\dagger}\hat{a}+\Delta\sum_{{\bm{i}}=1}^{L}{\hat{s}^{z}_{\bm{i}}}+2g_{a}(\hat{a}^{\dagger}+\hat{a})\sum_{\bm{i}=1}^{L}\hat{s}^{x}_{\bm{i}},$
(1)

where $L$ represents the number of TLSs in the system, $\Delta$ is the
difference in energy between ground and excited level in a two level systems
(TLS), and $\hat{a}^{\dagger}$ creates a photon in the cavity field.
Furthermore, the cavity mode frequency is $\omega_{a}$, $g_{a}$ represents the
coupling between a TLS and the cavity field (in the original paper [10],
$g_{a}\equiv\lambda/\sqrt{L}$, where $\lambda$ is used to characterise the
strength of the electron-photon coupling),
$\hat{s}^{j}_{\bm{i}}=\frac{1}{2}\sum_{\tau\tau^{\prime}}\hat{c}_{\tau,\bm{i}}^{\dagger}\sigma^{j}_{\tau\tau^{\prime}}\hat{c}_{\tau^{\prime},\bm{i}}$
are the spin operators, $\sigma^{j}$ are the Pauli spin matrices ($j=x,y,z$),
and $\hat{c}_{\tau,\bm{i}}^{\dagger}$, $\hat{c}_{\tau,\bm{i}}$ are the
electron creation and annihilation operators in $\bm{i}$th TLS. Index $\tau$
takes values $g$ (ground state) or $e$ (excited state).

Since we aim to investigate the SHG response from the Dicke model, the
Hamiltonian in Eq. (1) will be augmented with a second photon mode
representing a fluorescent field. Several studies devoted to the Dicke model
address the homogeneous case, where all the TLS are identical. There has
however been also considerable interest in studying disordered Dicke systems,
as they are relevant in describing situations such as cold atoms in a cavity
[16, 17]. In this respect, another way to enrich the physics of the model is
to consider the interaction between th electrons to make the system more
interesting.

The inclusion of interactions in the Dicke model has received some attention
in the literature, for example in the context of quantum information via the
so-called Hubbard-Dicke model for interacting qubits in coupled cavities (see
e.g. [18]), or interacting Rydberg atoms is optical cavities (see e.g. [19]).
In situations of this kind, mostly spin-spin type (for dipole-dipole, qubit-
qubit) interactions are considered between the relevant (atom-atom, qubit-
qubit, etc) units. This choice is not necessarily always physically motivated,
but in any case it makes these extended models amenable to simplified
numerical treatments.

The effect of both disorder and interactions on SHG in the Dicke model is the
target of the present investigation and we will consider a more general type
of interactions, namely Hubbard-type interactions between electrons in nearest
neighbour TLSs.

As method of choice to tackle large samples with disorder and this type of
interactions, we advocate the use of the non-equilibrium Green’s functions
(NEGF) method [20, 21, 22, 23, 24] within the so-called generalized Kadanoff-
Baym ansatz (GKBA) [25], and in a formulation apt to deal with electron-boson
systems [26, 27, 28]. Furthermore, to have a most transparent platform to
assess the scope of NEGF-GKBA, we will consider for simplicity only
interactions between electrons in the excited levels of the TLS. We will see
in the following that this model already represents a valuable benchmark for
NEGF-GKBA.

After introducing the model in detail in Sect. 2, we briefly review the NEGF-
GKBA method in Sect. 3, and then we proceed to a comparison of NEGF-GKBA
results with some exact diagonalization ED benchmarks in Sect. 4. This is
followed by a discussion of interaction and disorder effects on SHG in Sect.
5, and, in in Sect. 6, by an analysis of how cavity leakage from an imperfect
cavity manifests in the SHG spectra. Finally, after a brief discussion of the
third-harmonic generation response in Sect. 7, we conclude in Sect. 8 with
some conclusive remarks and a succinct outlook.

##  2 Model Hamiltonian

For convenience, we rewrite $\hat{H}_{DM}$ of Eq. (1) in an explicit electron
picture (in the rewriting, an obvious constant energy shift of the Hamiltonian
is neglected):

$\displaystyle\hat{H}_{DM}\rightarrow\hat{\tilde{H}}_{DM}=\epsilon_{g}\sum_{\bm{i}=1}^{L}\hat{c}_{g,\bm{i}}^{\dagger}\hat{c}_{g,\bm{i}}+\epsilon_{e}\sum_{\bm{i}=1}^{L}\hat{c}_{e,\bm{i}}^{\dagger}\hat{c}_{e,\bm{i}}+\omega_{a}\hat{a}^{\dagger}\hat{a}+g_{a}(\hat{a}^{\dagger}+\hat{a})\sum_{\bm{i}=1}^{L}(\hat{c}_{e,\bm{i}}^{\dagger}\hat{c}_{g,\bm{i}}+\hat{c}_{g,\bm{i}}^{\dagger}\hat{c}_{e,\bm{i}}),$
(2)

where $\epsilon_{g}\,(\epsilon_{e})$ corresponds to the ground (excited) level
energy in a TLS, $\hat{c}_{g/e,\bm{i}}$ destroys an electron in the
ground/excited level in $\bm{i}$th TLS (at all times, there is only one
electron in each TLS, and $\epsilon_{e}-\epsilon_{g}=\Delta$).

In this work we intend to explore the second harmonic generation (SHG) in a
Dicke system in the presence of e-e interactions and energies of the TLSs are
inhomogeneously distributed. To this end, $\hat{\tilde{H}}_{DM}$ is augmented
with a second photon mode, accounting for the fluorescent field. Furthermore,
to take into account e-e interactions, we consider the TLSs as arranged and
numbered according to a given sequence, and we use a notion of “adjacency”
between TLS with consecutive numbers. With these extensions and conventions,
the system’s Hamiltonian $\hat{\tilde{H}}_{DM}\rightarrow\hat{H}_{S}$, with

$\hat{H}_{S}=\hat{H}_{ele}+\hat{H}_{ph}+\hat{H}_{e-ph}.$ (3)

We now specify in detail the different terms in $\hat{H}_{S}$, starting from
the electronic part. This reads

$\displaystyle\hat{H}_{ele}=\sum_{\bm{i}}^{L}\tilde{\epsilon}_{g,\bm{i}}\hat{c}_{g,\bm{i}}^{\dagger}\hat{c}_{g,\bm{i}}+\sum_{\bm{i}}^{L}\tilde{\epsilon}_{e,\bm{i}}\hat{c}_{e,\bm{i}}^{\dagger}\hat{c}_{e,\bm{i}}+\sum_{<\bm{i},\bm{j}>}^{L}\frac{U_{e}}{2}\hat{n}_{e,\bm{i}}\hat{n}_{e,\bm{j}},$
(4)

The electronic interactions are introduced between excited electrons of TLSs
with nearest neighbour sites indexes ${\bm{i},\bm{j}}$ (denoted by the
$\langle{\bm{i},\bm{j}}\rangle$), with
$\hat{n}_{e,\bm{i}}=\hat{c}_{e,\bm{i}}^{\dagger}\hat{c}_{e,\bm{i}}$ the
density operator of the excited state of the ${\bm{i}}$-th TLS, and $U_{e}$
the strength of the interaction. A similar contribution could be considered
for the ground levels as well, but is omitted here for simplicity, and to
focus on the role of interactions in the system’s excited states. The effect
of inhomogeneity in the TLSs is included by choosing the ground and the
excited level energies as

$\displaystyle\tilde{\epsilon}_{{g,\bm{i}}}=\epsilon_{g}-\frac{\delta}{2}\sin\Big{[}\frac{\pi({\bm{i}}-1)}{L}\Big{]}$
(5)
$\displaystyle\tilde{\epsilon}_{{e,\bm{i}}}=\epsilon_{e}+\frac{\delta}{2}\sin\Big{[}\frac{3\pi({\bm{i}}-1)}{L}\Big{]}.$
(6)

These energy distributions will introduce a pseudo-disorder (in the following,
just referred to as “disorder” for simplicity) in the resonant frequencies of
the Dicke model, meaning that many TLSs are actually off-resonance whenever
$\delta\neq 0$, at an extent determined by the value of $\delta$.

The free term for the two photon fields are

$\hat{H}_{ph}=\omega_{a}\hat{a}^{\dagger}\hat{a}+\omega\hat{b}^{\dagger}\hat{b},$
(7)

where $\hat{b}$ destroys a fluorescent photon, whilst the interaction between
the electrons and the photons is represented by

$\hat{H}_{e-ph}=\left[g_{a}(\hat{a}^{\dagger}+\hat{a})+g_{b}(t)(\hat{b}^{\dagger}+\hat{b})\right]\sum_{{\bm{i}}}^{L}(\hat{c}_{{g,\bm{i}}}^{\dagger}\hat{c}_{{e,\bm{i}}}+\hat{c}_{e,\bm{i}}^{\dagger}\hat{c}_{g,\bm{i}}),$
(8)

where $g_{b}(t)=g^{\prime}e^{-\Gamma t}$. Here, $g^{\prime}$ is the initial
strength of the coupling between the electron and the fluorescent field, and
$\Gamma$ accounts via a phenomenological damping (broadening) effects left out
in the model, e.g., non-radiative transitions and/or mode leakages in the
cavity [7, 9, 29]. In what follows, we introduce the generalised indices
$(g,\bm{j})\rightarrow 2j-1$ and $(e,\bm{j})\rightarrow 2j$.

##  3 The NEGF method

A straightforward method to describe the behavior of our modified Dicke
Hamiltonian is exact diagonalization (ED), which is however limited to rather
small samples. Nevertheless, if one adopts a spin formalism for the Dicke
model (i.e. each TLS is represented as a spin 1/2 degree of freedom, then in
the absence of interaction and disorder the Hamiltonian commutes with the
total spin of the system. In that case, starting from a state of definite
total spin, the exact numerical time evolution can be performed by considering
only a spin-conserved subspace of the complete basis set. With this reduction,
the capacity of ED can be stretched to study larger Dicke systems. On the
other hand, with interaction and disorder present (in the sense specified in
Eq. 3) the total spin is not conserved, and ED can still be employed in the
electronic-level picture, but not for large samples.

This brings us to the central point of the paper, namely the proposal that for
large samples with disorder and/or interactions, we can make use of the NEGF-
GKBA in a formulation suitable for electron-boson systems [26, 27, 28].

Recently, a formulation of the NEGF-GKBA method was introduced [30, 31], which
attains time linear scaling in the simulation time. The approach has
subsequently been extended to case of interacting electron-boson systems [26,
27, 28]. In this work, we consider this linear-time formulation to study Dicke
systems of large size and in the rest of this section we provide a short
summary of the method, referring to the original literature for details.

In the theoretical treatment, the photon fields are represented by position
and momentum operators instead of creation/annihilation operators. For the
cavity field, the position and momentum operators are
$\hat{x}_{1}=\nobreak(\hat{a}^{\dagger}+\nobreak\hat{a})/\sqrt{2}$ and
$\hat{p}_{1}=i(\hat{a}^{\dagger}-\nobreak\hat{a})/\sqrt{2}$. Similarly, the
fluorescent field is represented by $\hat{x}_{2}$ and $\hat{p}_{2}$.
Furthermore, it is useful to represent jointly the position and momentum
operators via an operator $\hat{\phi}_{\mu}$, where the index $\mu$ is a
combined index with $\hat{\phi}_{i,1}=\hat{x}_{i}$ and
$\hat{\phi}_{i,2}=\hat{p}_{i}$. Hence, in terms of the
$\\{\hat{\phi}_{\mu}\\}$ operators, the Hamiltonian of the photon fields
results in
$\hat{H}_{ph}=\sum_{\mu\nu}\Omega_{\mu\nu}\hat{\phi}_{\mu}\hat{\phi}_{\nu}$
and, for the electron-photon interaction Hamiltonian, we can write, in the
same notation, $\hat{H}_{el-
ph}=\sum_{\mu,ij}g_{\mu,ij}(t)\hat{c}_{i}^{\dagger}\hat{c}_{j}\hat{\phi}_{\mu}$.

In the NEGF-GKBA scheme for electron-boson systems, the central ingredient is
the pair of density matrices $(\rho,{\bm{\gamma}})$, respectively for
electrons and for bosons [26, 27, 28]. More in detail, the (single-particle)
density matrix for electrons is
$\rho_{ij}(t)=\big{<}\hat{c}_{j}^{\dagger}(t)\hat{c}_{i}(t)\big{>}$, while the
boson density matrix is given by
$\gamma_{\mu\nu}(t)=\big{<}\Delta\hat{\phi}_{\nu}(t)\Delta\hat{\phi}_{\mu}(t)\big{>}$,
where $\Delta\hat{\phi}_{\nu}(t)\equiv\hat{\phi}_{\nu}(t)-\phi_{\nu}(t)$ and
$\phi_{\nu}(t)\equiv\langle\hat{\phi}_{\nu}(t)\rangle$. In terms of $\rho$ and
${\bm{\gamma}}$ it is then convenient to consider an effective one-particle
electronic Hamiltonian $h^{e}$ at the Hartree-Fock (HF) level,

$\displaystyle
h^{e}_{ij}(t)=h^{0}_{ij}+\sum_{m}\delta_{ij}\mathcal{U}_{im}\rho_{mm}^{<}(t)-\mathcal{U}_{ij}\rho_{ij}^{<}(t)+\sum_{\mu}\phi_{\mu}(t)g_{\mu,ij}(t),$
(9)

where $h^{0}_{ij}\hat{c}_{i}^{\dagger}\hat{c}_{j}$ denotes all the non
interacting contributions in Eq. (4), and the interaction terms are chosen as
$\mathcal{U}_{ij}=U_{ij}\delta_{i,j\pm 1}$ (the appropriate choice of $U_{ij}$
permits to specialise the type of interactions, e.g. only between excited
levels). The GKBA for fermions has been scrutinised quite in detail in the
literature. On the other hand, in this work we are primarily interested on the
scope of NEGF-GKBA for SHG in the multi-photon regime and for models of
quantum optics. This, together with our wish of avoiding unnecessarily heavy
calculations, motivates our simple, HF treatment of e-e interactions. In the
same spirit, the effective boson Hamiltonian $h^{b}_{\mu\nu}=\nobreak
2\sum_{\xi}\alpha_{\mu\xi}\Omega_{\xi\nu}$ with
$\alpha_{\mu\nu}=[\hat{\phi}_{\mu},\hat{\phi}_{\nu}]$.

The resulting equations of motion for the density matrices assume then the
form [27, 28]

$\displaystyle\frac{d}{dt}\rho^{<}(t)$
$\displaystyle=-i\big{[}h^{e}(t),\rho^{<}(t)\big{]}-(I^{e}(t)+{I^{e}}^{\dagger}(t)),$
(10) $\displaystyle\frac{d}{dt}{\bm{\gamma}}^{<}(t)$
$\displaystyle=-i\big{[}{\bm{h}}^{b}(t),{\bm{\gamma}}^{<}(t)\big{]}+({\bm{I}}^{b}(t)+{{\bm{I}}^{b}}^{\dagger}(t))$
(11)

where the collision integrals $I^{e}$ and ${\bm{I}}^{b}$ account for the
interaction between electrons and photon fields. More explicitly,
$I^{e}_{mj}=i\sum_{\mu,l}g_{\mu,ml}\mathcal{G}^{b}_{\mu,lj}$ and
$I^{b}_{\mu\nu}=-i\sum_{\eta,mj}\alpha_{\mu\eta}g_{\eta,mj}\mathcal{G}^{b}_{\nu,jm}$,
where the high-order Green’s function
$\mathcal{G}^{b}_{\mu,ij}(t)=\big{<}\hat{c}^{\dagger}_{j}(t)\hat{c}_{i}(t)\hat{\phi}_{\mu}(t)\big{>}_{c}$.

To explicitly compute the collision integrals within many-body perturbation
theory, use is then made of the electron and the boson GKBA scheme:

$\displaystyle G^{\lessgtr}(t,t^{\prime})$
$\displaystyle=-G^{R}(t,t^{\prime})\rho^{\lessgtr}(t^{\prime})+\rho^{\lessgtr}(t)G^{A}(t,t^{\prime})$
(12) $\displaystyle{\bm{D}}^{\lessgtr}(t,t^{\prime})$
$\displaystyle={\bm{D}}^{R}(t,t^{\prime}){\bm{\alpha}\bm{\gamma}}^{\lessgtr}(t^{\prime})-{\bm{\gamma}}^{\lessgtr}(t){\bm{\alpha}}{\bm{D}}^{A}(t,t^{\prime}),$
(13)

where $G$ denotes the single particle electron Green’s function, ${\bm{D}}$ is
the corresponding boson counterpart, $R$ ($A$) denote the retarded (advanced)
components, and the connection to the density matrices is given by
$G^{\lessgtr}(t,t)=i\rho^{\lessgtr}(t)$ and
${\bm{D}}^{\lessgtr}(t,t)=-i{\bm{\gamma}}^{\lessgtr}(t)$. In terms of the GKBA
expression above, a time linear scheme is obtained, where the ensuing coupled
ordinary differential equations read:

$\displaystyle i\frac{d}{dt}\phi_{\mu}$
$\displaystyle=\sum_{\nu}h^{b}_{\mu\nu}\phi_{\nu}+\sum_{\nu,ij}\alpha_{\mu\nu}g_{\nu,ij}\rho^{<}_{ji}$
(14a) $\displaystyle i\frac{d}{dt}\rho_{lj}^{<}$
$\displaystyle=[h^{e},\rho^{<}]_{lj}+\sum_{\mu,i}g_{\mu,li}\mathcal{G}_{\mu,ij}^{b}-(l\leftrightarrow
j)^{*},$ (14b) $\displaystyle i\frac{d}{dt}\gamma_{\mu\nu}^{<}$
$\displaystyle=[{\bm{h}}^{b},{\bm{\gamma}}^{<}]_{\mu\nu}+\sum_{\eta,mn}\alpha_{\mu,\eta}g_{\eta,mn}\mathcal{G}^{b}_{\nu,nm}-(\mu\leftrightarrow\nu)^{*},$
(14c) $\displaystyle i\frac{d}{dt}{\mathbfcal{G}}^{b}$
$\displaystyle=-{\bm{\Psi}}^{b}+{\bm{h}}^{b}{\mathbfcal{G}}^{b}-{\mathbfcal{G}}^{b}{\bm{h}}^{e},$
(14d)

where ${\bm{h}^{e}=h^{e}\otimes I-I\otimes(h^{e})^{T}}$,
${\bm{\Psi}}^{b}\equiv{\bm{\gamma}}^{>}{\bm{g}}{\bm{\rho}}^{<}-{\bm{\gamma}}^{<}{\bm{g}}{\bm{\rho}}^{>}$
with ${\bm{\rho}}^{>}=\rho^{>}\otimes(\rho^{<})^{T}$ and
${\bm{\rho}}^{<}=\rho^{<}\otimes(\rho^{>})^{T}$.

Within this NEGF-GKBA treatment, we use the average number of emitted photons
$\big{<}\hat{b}^{\dagger}\hat{b}\big{>}$ as indicator of the occurrence of SHG
and, as the system’s initial state, we take a product state of the electron
system and the two photon fields. As in e.g.[7, 9, 32, 29], initially the
electron system will be in the ground state, whilst the initial state of the
cavity field is a coherent state

$\displaystyle\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left|\beta\right\rangle$
$\displaystyle\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=e^{\beta\hat{a}^{\dagger}-\beta^{*}\hat{a}}|0\rangle,$
(15)

and the emitted photon field is initially represented by the vacuum state
$\left|0\right\rangle$. In all the calculations, the resonance frequency of a
TLS without disorder is taken ${\Delta}=\epsilon_{e}-\epsilon_{g}=1.0$, the
cavity field coupling is $g_{a}=0.03$, and the fluorescent field coupling
$g^{\prime}=0.01$, with $\Gamma=0.02$. The coherent state of the cavity field
is chosen to contain
$n_{a}=\langle\hat{a}^{\dagger}\hat{a}\rangle=\beta^{2}=9$ average number of
photons. Finally, the cavity frequency in all the SHG calculations is
${\omega_{a}}={\Delta}/2=0.5$.

##  4 ED and NEGF-GKBA

To benchmark the NEGF-GKBA approach, we have considered ED calculations for a
smaller number of TLS, and compared these with GKBA results. In most cases,
both disorder (in the sense specified in Eq. (5)) and e-e interactions are
included. The comparisons are shown in Figure 1, where it is manifest that
overall ED and GKBA results agree very well. Good agreement is observed also
in the presence of disorder. However, with e-e interactions included, the SHG
peak from GKBA is shifted slightly towards the higher frequencies compared to
ED.

$\bm{(a)}$$\bm{(b)}$$\bm{(c)}$$\bm{(d)}$$\bm{(e)}$$\bm{(f)}$$L=3$$L=7$ Figure
1: SHG spectrum from ED (brown curves) and NEGF-GKBA (solid filled blue
curves) for $L=3$ (upper panel) and $L=7$ (lower panel). The results are for
disorder strengths $\delta=0.0\,(\bm{a},\bm{d}),\,\delta=0.1\,(\bm{b},\bm{e})$
and $\delta=0.2\,(\bm{c},\bm{f})$, and interaction $U_{e}=0.0,0.3$ and $0.5$.
All curves are shown at $t_{e}=250$, the end point of the time-evolution
interval.

To characterize these discrepancies in more detail, we singled out the case of
three TLS (i.e., $L=3$) for several values of $U_{e}$ and, to establish
trends, we focused on the behavior of the maximum intensity in the SHG peak
$\langle\hat{b}^{\dagger}\hat{b}\rangle_{max}$ in relation of the
corresponding emitted frequency $\omega_{max}$. This was done for both ED and
GKBA, and the inherent results are shown in Figure 2$\bm{(a)}$ and
2$\bm{(b)}$. In ED, at larger $U_{e}$, $\omega_{max}$ increases only slightly
(panel $\bm{(a)}$), i.e. the SHG peak appears to be not very sensitive to
$U_{e}$ (at least in this interaction range). However, with GKBA, a larger
shift of the SHG peak occurs with the increase of $U_{e}$.

Coming now to the SHG spectral intensity, in panel $\bm{(b)}$ we observe that,
in the ED calculations, increasing $U_{e}$ reduces the intensity of the SHG
peak, whereas, in GKBA, the trend of decreasing intensity is observed only for
lower $U_{e}$; indeed, for higher $U_{e}$ the intensity almost remains
constant. It is beyond the scope of this short contribution to address in a
quantitative way the reasons for these small but clearly discernible
discrepancies. Still, it can be of interest to provide some additional
considerations. A first factor to take into account is that, together with
photon dressing effects, the system’s energy levels undergo multiple level
crossing as a function of the interaction $U_{e}$ ( in particular, the level
crossing at $U_{e}=0$ represents an essential complication, which is not easy
to account in perturbation theory. The second aspect to consider is that, in
any case, an HF treatment of e-e interactions can be inadequate, as shown next
with simple qualitative argument in the large-$U_{e}$ limit.

Specifically, in our model (where interactions only occur among electrons in
the upper levels of the TLSs), the ground state within an HF description is a
state with, say, a fraction $\delta$ of electrons on the upper level. This
fraction can be estimated by balancing the energy between the electrons and
the $\omega_{a}$-type photons:
$U_{e}\delta^{2}\approx{\omega_{a}\beta^{2}\rightarrow\delta\approx\sqrt{\omega_{a}\beta^{2}/U_{e}}}$.
Excited states in the HF theory are rather simple and similar to the exact
ones (in our case we basically have a spin-flip transition of an electron from
the lower to the upper state in the TLS). The energy of such a transition, in
addition to the noninteracting part equal to 1, contains a mean field
contribution
$U_{e}[(1-\delta)^{2}-\delta^{2}]=U_{e}(1-2\delta)={U_{e}[1-2\sqrt{\omega_{a}\beta^{2}/U_{e}}]}$.
This yields in total for the fluorescence $\omega\approx
1+U_{e}[1-2\sqrt{\omega_{a}\beta^{2}/U_{e}}]$, in contrast to $\omega\simeq 1$
in the exact theory.

Even so, it is important to restate here that the overall performance of the
NEGF-GKBA is very good, together with the observation (stimulated by
preliminary estimates, not shown here) that the differences observed most
likely are due to the shortcomings of the HF treatment for the e-e
interactions, and not of the photon-electron interactions.

$\bm{(a)}$$\bm{(b)}$${\bm{(c)}}$ Figure 2: ED and NEGF-GKBA results at $L=3$
and $\delta=0.0$ for the frequency $\omega_{max}$ (panel $\bm{(a)}$) and
photon occupation $\langle\hat{b}^{\dagger}\hat{b}\rangle$ (panel $\bm{(b)}$)
in correspondence of the maximum intensity of the SHG spectrum. Panel
${\bm{(c)}}$: ED (black empty curve) vs GKBA (solid filled curves) results for
Dicke systems at $\delta=U_{e}=0.0$ and $L=10,20,30,40,50$ and $60$.In all
panels, the ED calculations are performed with the spin formalism, and all
curves are calculated at time $t_{e}=250$.

###  4.1 Spin representation in ED for large Dicke systems

The comparisons performed so far, while successful for NEGF-GKBA, concerned
rather small systems. However, a main aim of this work is to explore the scope
of NEGF-GKBA for large and interacting and/or disordered Dicke systems, in
general not accessible to ED. Yet, at least in the absence of disorder and e-e
interactions, it is still possible to use ED as reference for $L$ considerably
larger than in the previous section. To make this explicit, we consider again
the spin operators representation so that (dropping again the irrelevant total
energy shift in  $\hat{H}_{S}$) the Hamiltonian in Eq. (3) becomes

$\hat{\tilde{H}}_{S}={\Delta}\hat{S}^{z}+\omega_{a}\hat{a}^{\dagger}\hat{a}+\omega\hat{b}^{\dagger}\hat{b}+2\left[g_{a}(\hat{a}^{\dagger}+\hat{a})+g_{b}(t)(\hat{b}^{\dagger}+\hat{b})\right]\hat{S}^{x},$
(16)

where $\hat{S}^{x,z}=\sum_{\bm{i}=1}^{L}\hat{s}_{\bm{i}}^{x,z}$.

In this form, $\hat{\tilde{H}}_{S}$ manifestly commutes with the total spin
$\hat{S}^{2}$ operator, and thus the photon-induced transitions conserve the
total spin of the system. Hence, starting from an initial state of definite
total spin (e.g., with the electron in each TLS in the lower level,
corresponding to $S=L/2$ and $S_{z}=-L/2$), the dynamics of the non-disordered
and non-interacting electron system induced by the photon fields can be
represented in ED within a Hilbert space of $L+1$ (instead of $2^{L}$) basis
states. This makes possible to consider ED benchmarks (for $U_{e}=\delta=0$)
when $L$ is large, and thus have a further comparison with NEGF-GKBA. The
comparison is shown in Figure 2$\bm{(c)}$, with a black empty curve
representing ED results, and the coloured solid curve corresponds to GKBA
results. We can observe that ED and GKBA results are in excellent agreement in
the considered parameter range. Commenting on the physical information in
these results, we note that until $L=40$, the intensity of the spectra
increases with $L$, which is expected. However, for $L>40$, even though there
is a small increase in the intensity of SHG with $L$, the intensity of the
main peak (peak close to $\omega=0.5$) decreases with $L$.

Figure 2$\bm{(c)}$ shows the time evolved plots obtained at a final time
$t_{e}=250$. It can be useful to consider the actual development in time of
such spectra, that is reported in Figure 3. During the time evolution, and for
$L<20$, the intensity of the spectra increases initially but then remains
almost constant. However, for larger Dicke systems ($L>20$), the trend is not
monotonic: even though we initially observe a rise in the intensity in time,
at some point the signal will start to drop, which is consistent with the
results of Figure 2$\bm{(c)}$ in the long time limit. In summary, for the
parameters considered, a steady increase in the intensity of the SHG spectra
on increasing $L$ was not observed.

$\bm{(a)}$$\bm{(b)}$$\bm{(c)}$$\bm{(d)}$ Figure 3:  Time dependent NEGF-GKBA
spectra for $\delta=U_{e}=0.0$ and $L=10,20,30,60$. The parameters used are
$\epsilon_{e}-\epsilon_{g}\equiv\Delta=1.0,g_{a}=0.03,g^{\prime}=0.01,\hat{n}_{a}=\langle\hat{a}^{\dagger}\hat{a}\rangle=\beta^{2}=9,\Gamma=0.02$,
and $\omega_{a}=0.5$ All curves are shown at $t_{e}$=250, the end point of the
time-evolution interval.

##  5 Disorder and interaction

We now consider the role of disorder and electron interactions on SHG for the
larger Dicke systems, and only within a NEGF-GKBA treatment. Disorder and
interactions are included either one at the time (Figure 4), to disentangle
their effect, or together (Figure 5), to see if they are mutually competing or
synergic in affecting the spectrum.

SHG spectra for $U_{e}=0$ and different disorder strengths (Figure
4$\bm{(a)}$) show that at larger $\delta$ the SHG spectral will diminish,
along with some peak splitting. Increasing $\delta$ will set the resonance
frequency ($\omega_{res}$) of several TLS: in the Dicke system away from
$\Delta$. Hence the SHG peak intensity decreases. Also, the different TLS of
the disordered system will have different resonance frequencies; hence, there
are multiple peaks around the nominal SHG frequency. Conversely, in Figure
4$\bm{(b)}$ there is no disorder ($\delta=0.0$), and one can observe the sole
effect of electron interaction. On increasing the electron interaction, the
intensity of the SHG peak reduces again. But the intensity reduction observed
with the electron interaction is less compared to the disorder induced
reduction. Finally, we also considered calculations with both disorder and
interactions present (Figure 5). We did not observe any mutual
compensation/competition of disorder and interactions effects; rather, the
intensity of the SHG peak got reduced even further.

$\bm{(a)}$$\bm{(b)}$ Figure 4: Panel $\bm{(a)}$: spectra in the presence of
disorder, with $\delta=0.0,0.1,0.2$ and $U_{e}=0.0$. Panel $\bm{(b)}$: spectra
in the presence of interactions, with $U_{e}=0.0,0.3,0.5$ and $\delta=0.0$.
The systems considered have $L=10,15,20,25,30$. The spectra for $L=10$ are
scaled by ‘$3$’ and $L=15$ by ‘$2$’ for visual clarity. All spectra were
obtained at the time $t_{e}=250$. $\bm{(a)}$$\bm{(b)}$$\bm{(c)}$$\bm{(d)}$
Figure 5:  Fluorescent spectra obtained with NEGF-GKBA for $L=10,15,20,25$,
disorder strength $\delta=0.0$ and interaction strength $U_{e}=0.0$. The
spectra were calculated at the time $t_{e}=250$.

##  6 Cavity leakage

In reality, photons inside the cavity will escape from the cavity after some
time. In the case of a bad cavity, the photons will escape more quickly. We
now mimic the case of a bad cavity by introducing cavity leakage via coupling
the photon fields to a bath of classical harmonic oscillators. This method is
inspired from the physics of the Caldeira–Leggett model [33, 34]. Hence, the
system Hamiltonian $\hat{H}_{S}$ in Eq. 3 will be augmented with the two extra
terms $\hat{H}_{ph-bath}$ and $\hat{H}_{bath}$ [29]. The interaction between
the photon field and the bath is represented by

$\hat{H}_{ph-
bath}=\sum_{k}^{N_{bath}}C_{k}\left[{q_{a,k}}(t)\big{(}\hat{a}^{\dagger}+\hat{a}\big{)}+{q}_{b,k}(t)\big{(}\hat{b}^{\dagger}+\hat{b}\big{)}\right],$
(17)

where ${q_{a,k}}$ and ${{q}_{b,k}}$ are the coordinates of the oscillators,
the coupling strength is $C_{k}$ and $C_{k}=A({\Delta_{B}}k)^{a}$.

The Hamiltonian  function corresponding to the baths of the classical
oscillators is

$H_{bath}=\sum_{k}^{N_{bath}}\Big{(}\frac{{p_{a,k}^{2}}+{p_{b,k}^{2}}}{2m_{k}}+\frac{1}{2}\omega_{k}^{2}m_{k}\big{[}{q_{a,k}^{2}}+{{q}_{b,k}^{2}}\big{]}\Big{)},$
(18)

where $p_{a,k}\,(p_{b,k})$ represents the momentum of the $k$th oscillator
coupled to the cavity (fluorescent) field, and $\omega_{k}={\Delta_{B}}k$
corresponds to the frequency of the oscillator.

In terms of the operator $\hat{\phi}_{\mu}$ defined in Sect. 3, one can
express the quantum-classical Hamiltonian as $\hat{H}_{ph-
bath}=\sum_{\mu}F_{\mu}(t)\hat{\phi}_{\mu}$ where, explicitly,
$\\{F_{\mu}\\}\equiv(F_{1,1},F_{1,2},F_{2,1},F_{2,2})$, with
$F_{1,1}=\sqrt{2}\sum_{k}^{N_{bath}}C_{k}{q_{a,k}}(t)$,
$F_{2,1}=\sqrt{2}\sum_{k}^{N_{bath}}C_{k}{q_{b,k}}(t)$ and
$F_{1,2}=F_{2,2}=0$.

In the presence of cavity leakage, the NEGF-GKBA scheme in the Eqs. (14a-d)
will remain the same except for Eq. 14(a) [24], since now the time evolution
of the average position will have a contribution from the cavity leakage:

$i\frac{d}{dt}\phi_{\mu}=\sum_{\nu}h^{b}_{\mu\nu}\phi_{\nu}+\sum_{\nu,ij}\alpha_{\mu\nu}g_{\nu,ij}\rho^{<}_{ji}+\sum_{\nu}\alpha_{\mu\nu}F_{\nu}.$
(19)

Thus, together with the NEGF-GKBA time-evolution, the coordinates of the
classical oscillators are time-evolved using the Ehrenfest dynamics [24, 35].
In the calculations, we considered all the oscillators’ masses $m_{k}=1$.
Hence,

$\displaystyle{\ddot{q}_{a,k}}(t)$
$\displaystyle=\nobreak-\omega_{k}^{2}{q_{a,k}}(t)+C_{k}\langle\hat{a}^{\dagger}+\hat{a}\rangle_{{\bar{q}_{a},\bar{q}_{b},t}},$
(20) $\displaystyle{\ddot{{q}}_{b,k}}(t)$
$\displaystyle=\nobreak-\omega_{k}^{2}{{q}_{b,k}}(t)+C_{k}\langle\hat{b}^{\dagger}+\hat{b}\rangle_{{\bar{q}_{a},\bar{q}_{b},t}},$
(21)

where $\bar{q}_{a}=\\{{q_{a,k}}\\}$ and $\bar{q}_{b}=\\{{q_{b,k}}\\}$, also,
we have considered that, at initial time, the bath oscillators are in their
equilibrium positions, i. e. $q_{a,k}=q_{b,k}=0$ and $p_{a,k}=p_{b,k}=0$. The
effects of cavity leakage are shown in Figure 6, where the pink curve
corresponds to the long-time limit result without cavity leakage, and the
time-evolving plots are obtained by including cavity leakage (the values of
the parameters used in the calculations are reported in the figure caption).
In all the cases shown, the time-evolved plots are less intense than the
corresponding pink curve. That is, cavity leakage reduces the intensity of the
fluorescent spectra also in the presence of disorder and interaction. This is
expected, since leakage removes photons from the optical cavity.

$\bm{(a)}$$\bm{(b)}$$\bm{(c)}$$\bm{(d)}$$\bm{(e)}$$\bm{(f)}$$L=15$$L=25$
Figure 6: The effect of cavity leakage. The time-dependent plots (the pink
curves) refer to results with (without) the classical baths for $L=15$ (upper
panels) and $L=25$ (lower panels). The values of $\delta$ and $U_{e}$ used in
a given panel are indicated at the top of the panel. The bath parameters for
the cavity leakage calculations are $A=0.005$, $a=0.6$, ${\Delta_{B}}=0.01$
and $N_{bath}=200$. The time evolving plots are scaled by ‘$2$’ for visual
clarity.

##  7 Third harmonic generation

All the results so far were for the second harmonic spectrum, i.e. we
considered ${\omega_{a}}={\Delta}/2$. Clearly, the same NEGF-GKBA (and ED)
treatments can be employed to study third harmonic generation (THG), with
${\omega_{a}}={\Delta}/3$. In comparison to SHG, the THG peaks are less
intense, as it can be observed in Figure 7 $\bm{(a)}$, where the results were
obtained with ED. We then compared THG spectra obtained with the ED and NEGF-
GKBA in the case of non-interacting, non-disordered system. As shown in Figure
7$\bm{(b)}$, there is a quite good agreement between the results from the two
methods.

$\bm{(a)}$$\bm{(b)}$ Figure 7:  Panel $\bm{(a)}$: comparison of fluorescent
spectra for SHG and THG obtained with ED. Panel $\bm{(b)}$: comparison between
THG fluorescent spectra obtained with ED and NEGF-GKBA. All the SHG (THG)
calculations in figure are with ${\omega_{a}}={\Delta}/2=1/2$
(${\omega_{a}}={\Delta}/3=1/3$). For both panels, $U_{e}=\delta=0.0$, with the
spectra obtained at time $t_{e}=250$.

Finally, THG results from NEGF-GKBA (and their comparisons to ED results) when
disorder and interaction are included are shown in Figure 8. Again, the
agreement between the two approaches is quite satisfactory. Also, as for SHG,
both disorder and interactions will reduce the intensity of THG (the damping
of THG observed with the disorder present is larger than with the
interaction). As observed already for the case of SHG, the THG peak obtained
with NEGF-GKBA for $U_{e}=0.3$ is slightly shifted compared to the ED peak. We
expect the source of this discrepancy to be same as for SHG. In these very
preliminary discussion of THG, we examined only very small Dicke systems.
However, from what was found, and as suggested by the picture emerged for SHG,
it would appear that the NEGF-GKBA method can be successfully used to study
THG spectra in large Dicke systems, and in the presence of disorder and
interaction.

$\bm{(a)}$$\bm{(b)}$ Figure 8:  Comparison between THG fluorescent spectra
from ED and NEGF-GKBA treatments, and with or without disorder and electron
interactions (the used values of $\delta$ and $U_{e}$ are shown in the
legends).The results are for small systems ($L=3,4$), and the spectra are
obtained at time $t_{e}=250$.

##  8 Conclusions

The Dicke model is one of the paradigmatic system in quantum optics to address
nonlinear light-matter interaction quantum phenomena. In this paper, we mainly
studied the second-harmonic generation (SHG) fluorescent spectra of the Dicke
model in the presence of disorder and electron interaction. To describe large
disordered and interacting samples, we proposed the use of nonequilibrium
Green’s functions (NEGF) in the so-called generalized Kadanoff-Baym ansatz
(NEGF-GKBA) for electron-boson systems, a state-of-art approach based on
nonequilibrium many-body perturbation theory.

The scope of the approach to describe SHG in the Dicke model was tested by
comparing NEGF-GKBA and exact diagonalization (ED) spectra for small
interacting disordered systems and larger clean noninteracting ones. For the
parameter regimes investigated, we observed a very good agreement between the
ED and NEGF-GKBA treatments (though, on increasing the electron-electron
interactions, the SHG peaks are slightly shifted in NEGF-GKBA as compared to
ED). These findings point at NEGF-GKBA as a powerful, rather accurate and
viable method to deal with quantum systems of some complexity. Furthermore

Concerning specific physical trends, our study suggests that both disorder and
electron-electron interactions cause a reduction in SHG intensity in the Dicke
model, and that disorder hinders SHG more strongly than interaction. We also
considered the role of cavity leakage by coupling the photon fields to baths
of classical oscillators. The results showed (or, more properly speaking,
confirmed) that, in accordance with one’s intuition, cavity leakage depresses
the intensity of the fluorescent spectra.

Further, we briefly explored the application of the approach to the
investigation of third harmonic generation (THG). We considered only small
samples, but these endorsed quite clearly the use of the method for THG.
Again, for THG, results from ED and GKBA agree very well, and similarly to
SHG, the THG signal decreases with disorder and interactions. That is, the
method can be applied also to the study of higher harmonic generation
processes.

Overall, our outlook expectation (and in part speculation) is that the NEGF-
GKBA approach is generally robust to deal with non linear phenomena in the
Dicke model and/or related models of quantum optics (and their
generalisations), and that a mean-field treatment of e-e correlations is a
useful starting point. However, to deal with degenerate states and range of
parameters in the vicinity of energy level crossings, higher order
diagrammatic contributions are to be included in the NEGF-GKBA self-energies.

Funding information: This research performed in this work was part of projects
funded by the Swedish Research Council with grants 2017-03945 and 2022-04486
(M.G. and C.V.) This research was part of project no. 2021/43/P/ST3/03293
cofunded by the National Science Centre and the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement no. 945339 (Y.P.).

## References

  * [1] N. Bloembergen,  _Rev. Mod. Phys._ 1982, _54_ , 3 685.
  * [2] M. F. Ciappina, J. A. Pérez-Hernández, A. S. Landsman, W. A. Okell, S. Zherebtsov, B. Förg, J. Schötz, L. Seiffert, T. Fennel, T. Shaaran, et al.,  _Rep. Prog. Phys._ 2017, _80_ , 5 054401.
  * [3] X. Fu, T. J. Cui,  _Progress in Quantum Electronics_ 2019, _67_ 100223\.
  * [4] C. Andraud, O. Maury,  _Eur. J. Inorg. Chem._ 2009, _2009_ , 29-30 4357.
  * [5] S. Yue, M. N. Slipchenko, J.-X. Cheng,  _Laser & photonics reviews_ 2011, _5_ , 4 496.
  * [6] G. F. Combes, A.-M. Vučković, M. Perić Bakulić, R. Antoine, V. Bonačić-Koutecky, K. Trajković,  _Cancers_ 2021, _13_ , 16 4206.
  * [7] M. Cini, A. d’Andrea, C. Verdozzi,  _Phys. Lett. A_ 1993, _180_ , 6 430.
  * [8] M. Cini, A. d’Andrea, C. Verdozzi,  _Int. J. Mod. Phys. B_ 1995, _9_ , 10 1185.
  * [9] E. Viñas Boström, A. D’Andrea, M. Cini, C. Verdozzi,  _Phys. Rev. A_ 2020, _102_ , 1 013719.
  * [10] R. H. Dicke,  _Physical review_ 1954, _93_ , 1 99.
  * [11] P. Kirton, M. M. Roses, J. Keeling, E. G. Dalla Torre,  _Advanced Quantum Technologies_ 2019, _2_ , 1-2 1800043\.
  * [12] E. Jaynes, F. Cummings,  _Proceedings of the IEEE_ 1963, _51_ , 1 89.
  * [13] H. J. Carmichael, D. F. Walls,  _J. Phys. B: At. Mol. Opt. Phys_ 1976, _9_ , 8 1199\.
  * [14] M. Tavis, F. W. Cummings,  _Phys. Rev._ 1968, _170_ , 2 379.
  * [15] J. Keeling,  _Phys. Rev. A_ 2009, _79_ , 5 053825.
  * [16] P. R. Eastham, P. B. Littlewood,  _Phys. Rev. B_ 2001, _64_ , 23 235101.
  * [17] O. Tsyplyatyev, D. Loss,  _Phys. Rev. A_ 2009, _80_ , 2 023803.
  * [18] S. He, L.-W. Duan, Y.-Z. Wang, C. Wang, Q.-H. Chen,  _Phys. Rev. A_ 2022, _106_ , 3 033712.
  * [19] X.-F. Zhang, Q. Sun, Y.-C. Wen, W.-M. Liu, S. Eggert, A.-C. Ji, et al.,  _Phys. Rev. Lett._ 2013, _110_ , 9 090402.
  * [20] L. Kadanoff, G. Baym,  _Quantum statistical mechanics Green’s function methods in equilibrium and nonequilibrium problems_ ,  W.A. Benjamin, New York, 1962.
  * [21] L. V. Keldysh,  _Sov. Phys. JETP_ 1965, _20_ 1018\.
  * [22] K. Balzer, M. Bonitz,  _Nonequilibrium Green’s function approach to inhomogeneous systems_ ,  Number 867 in Lecture notes in physics. Springer, Heidelberg, 2013.
  * [23] G. Stefanucci, R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ ,  Cambridge University Press, Cambridge, 2013.
  * [24] N. Säkkinen, Y. Peng, H. Appel, R. van Leeuwen,  _J. Chem. Phys._ 2015, _143_ , 23 234102.
  * [25] P. Lipavský, V. Špička, B. Velický,  _Phys. Rev. B_ 1986, _34_ , 10 6933.
  * [26] D. Karlsson, R. van Leeuwen, Y. Pavlyukh, E. Perfetto, G. Stefanucci,  _Phys. Rev. Lett._ 2021, _127_ , 3 036402.
  * [27] Y. Pavlyukh, E. Perfetto, D. Karlsson, R. van Leeuwen, G. Stefanucci,  _Phys. Rev. B_ 2022, _105_ , 12 125134.
  * [28] Y. Pavlyukh, E. Perfetto, G. Stefanucci,  _Phys. Rev. B_ 2022, _106_ , 20 L201408.
  * [29] M. Gopalakrishna, E. Viñas Boström, C. Verdozzi,  _SciPost Physics_ 2023, _15_ , 4 138.
  * [30] N. Schlünzen, J.-P. Joost, M. Bonitz,  _Phys. Rev. Lett._ 2020, _124_ , 7 076601.
  * [31] J.-P. Joost, N. Schlünzen, M. Bonitz,  _Phys. Rev. B_ 2020, _101_ , 24 245101.
  * [32] Y. Pavlyukh, E. Perfetto, D. Karlsson, R. van Leeuwen, G. Stefanucci,  _Phys. Rev. B_ 2022, _105_ , 12 125135.
  * [33] A. O. Caldeira, A. J. Leggett,  _Annals of physics_ 1983, _149_ , 2 374.
  * [34] V. Venkataraman, A. D. K. Plato, T. Tufarelli, M. S. Kim,  _J. Phys. B: At. Mol. Opt. Phys_ 2013, _47_ , 1 015501\.
  * [35] E. Boström, M. Hopjan, A. Kartsev, C. Verdozzi, C.-O. Almbladh,  _J. Phys. Conf. Ser._ 2016, _696_ , 1 012007.