Present address: ]MRC London Institute of Medical Sciences, Du Cane road,
London, W12 0NN, United Kingdom Present address: ]Centre for Nanophotonics,
Department of Physics, Engineering Physics & Astronomy, 64 Bader Lane, Queen’s
University, Kingston, Ontario, Canada K7L 3N6

# Sub-radiant states for imperfect quantum emitters coupled by a nanophotonic
waveguide

Xiao-Liu Chu  [<EMAIL_ADDRESS>Vasiliki Angelopoulou  Center for Hybrid
Quantum Networks (Hy-Q), Niels Bohr Institute, University of Copenhagen,
Blegdamsvej 17, DK-2100 Copenhagen, Denmark
Peter Lodahl  Nir Rotenberg  [  Center for Hybrid Quantum Networks (Hy-Q),
Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, DK-2100
Copenhagen, Denmark

###### Abstract

Coherent interactions between quantum emitters in tailored photonic structures
is a fundamental building block for future quantum technologies, but remains
challenging to observe in complex solid-state environments, where the role of
decoherence must be considered. Here, we investigate the optical interaction
between two quantum emitters mediated by one-dimensional waveguides in a
realistic solid-state environment, focusing on the creation, population and
detection of a sub-radiant state, in the presence of dephasing. We show that
as dephasing increases, the signatures of sub-radiance quickly vanish in
intensity measurements yet remain pronounced in photon correlation
measurements, particularly when the two emitters are pumped separately so as
to populate the sub-radiant state efficiently. The applied Green’s tensor
approach is used to model a photonic crystal waveguide, including the
dependence on the spatial position of the integrated emitter. The work lays
out a route to the experimental realization of sub-radiant states in
nanophotonic waveguides containing solid-state emitters.

††preprint: APS/123-QED

##  I Introduction

With the maturation of quantum photonic platforms, light-matter interaction is
no longer limited to the control of single emitters [1, 2, 3], as multi-
emitter systems are increasingly being developed and investigated [4, 5, 6,
7]. Of the different platforms, waveguide quantum electrodynamics (w-QED)
where quantum emitters are coupled to photonic waveguides [8, 4, 9, 5, 6, 10],
is particularly promising. Solid-state w-QED systems allow on-chip integration
based on reliable nanofabrication [11, 9, 12], permitting electrical control
over the emitter resonances [13, 2, 14], although scaling-up to many emitters
is challenged by inhomogeneous broadening unlike the case of atoms [10].
Importantly, even few quantum emitters deterministically coupled to a
waveguide can be a very powerful quantum resource since each emitter can
produce a high number of photonic qubits. Indeed two coherently-coupled
quantum emitters have been proposed as sufficient for creating large-scale
photonic cluster states for quantum communication [15, 16] or to access
decoherence-free sub-spaces for quantum computation [17]. The latter requires
the creation of long-lived sub-radiant states: collective excitations with
lifetimes that are much longer than those of the single emitters.

Previous studies have focused on systems of ideal emitters [18, 19, 20], while
here we show how sub-radiant states may be created and populated using real-
world quantum emitters in the presence of imperfections. We specifically
consider solid-state emitters [21] coupled to nanophotonic waveguides and
focus on the role of pure dephasing. We show that while even a small amount
dephasing (less than $10\%$ of the emitter decay rate $\Gamma_{0}$) is
sufficient to remove the signature of sub-radiance in photonic intensity
measurements, the states can be efficiently excited and detected from photon
correlation measurements. We observe that dephasing introduces a mixing
between the super and sub-radiant states, providing another pathway to
populate states that would otherwise be inaccessible. Finally, we calculate
the Green’s tensor for a photonic crystal waveguide and apply that to
determine the spatial dependence of the emitter-emitter coupling.

Figure 1: a. Two spatially separated emitters (by $\Delta z$) couple to the
same electromagnetic field of a single waveguide mode (See inset). This
coupled system may be excited either through the waveguide mode (resonant
transmission) or from free-space (resonant fluorescence), in which case the
emitters may be addressed individually. The resultant four energy levels of
the coupled system when the coupling is b) fully dispersive or c) fully
dissipative, along with d) exemplary time-dependent emission traces from the
different transitions. For dispersive coupling, the decay rates of the coupled
transitions remain the same as that of the individual systems, corresponding
to the green curve in d), while the energy levels are shifted. In contrast,
for dissipative coupling the energy levels are unchanged while the decay rates
are altered by $\Gamma_{12}$ such that one transition is sub-radiant (shown in
red here and in d) and on is super-radiant (shown in blue here and in d).
Depending on emitter position, the coupling may be partially dispersive and
partially dissipative.

##  II  sub-radiance of two quantum emitters in a nanophotonic waveguide

###  II.1 System geometry, excitation and measurements

The overall goal is to exploit how to efficiently control and populate a sub-
radiant state using coupled quantum emitters in a nanophotonic waveguide and
study the robustness towards imperfections. The system is sketched in Fig. 1
and is comprised of two, two-level emitters (TLEs) with identical resonant
frequencies $\omega_{0}\equiv\omega_{1}=\omega_{2}$, and decay rates
$\Gamma_{0}\equiv\Gamma_{1}=\Gamma_{2}$. Both emitters couple to the same
guided mode, and are separated by a distance $\Delta
z_{12}\gg\lambda_{\mathrm{wg}}$. Here, $\lambda_{\mathrm{wg}}$ is the
effective wavelength of the guided mode, at $\omega_{0}$, meaning that the
TLEs are coupled via the long-range dipole-dipole interaction mediated by the
guided mode [19].

When coupled, excitations are no longer ascribed to individual emitters, but
rather the system as a whole, resulting in the Dicke energy states as shown in
Fig. 1b and c. These coupled states are [22]:

$\displaystyle\ket{g}=\ket{g_{1}}\otimes\ket{g_{2}},$
$\displaystyle\ket{s}=\frac{1}{\sqrt{2}}(\ket{e_{1}}\otimes\ket{g_{2}}+\ket{g_{1}}\otimes\ket{e_{2}}),$
$\displaystyle\ket{a}=\frac{1}{\sqrt{2}}(\ket{e_{1}}\otimes\ket{g_{2}}-\ket{g_{1}}\otimes\ket{e_{2}}),$
$\displaystyle\ket{e}=\ket{e_{1}}\otimes\ket{e_{2}}.$ (1)

Not only are the intermediate states $\ket{s}$ and $\ket{a}$ entangled [23,
24], but both the decay rates and energy levels are modified by
$\Gamma_{12}=\Gamma_{21}$ and $J_{12}=J_{21}$, respectively, as shown in the
energy diagram (the full expression for these terms is detailed in the next
section). In what follows, we focus on identifying and populating the sub-
radiant state, whose decay rate
$\Gamma_{\mathrm{s/a}}=\Gamma_{0}\pm\Gamma_{12}\ll\Gamma_{0}$. In the case
where $\Gamma_{12}$ is a positive quantity (negative), the asymmetric state
$\ket{a}$ is sub-radiant (and $\ket{s}$ super-radiant), and vice versa for a
negative $\Gamma_{12}$.

There are two ways to excite the w-QED system, both of which are shown in Fig.
1a. Either both emitters are excited through the waveguide in the resonant
transmission (RT) configuration, or they can be addressed through free-space
in the resonant fluorescence (RF) configuration. For RT, the relative phases
with which the emitters are excited is locked by the properties of the guided
mode to which they couple, and determined by their separation. In this
configuration, the transmitted signal includes pump photons that coherently
interfere with those scattered from the emitters [3]. In contrast, free-space
excitation allows either, or both, of the emitters to be addressed with any
arbitrary phase or intensity difference. Here, only photons that originate
from the emitters are detected.

Regardless of the excitation method, we detect the light field at the output
port of the waveguide, as shown in Fig. 1a. In practice this requires out-
coupling of the light, for example using a grating as shown in the figure,
while we calculate the field in the waveguide far away from the emitters. We
calculate the intensity spectrum in either of the two configurations,

$\displaystyle\text{I}_{\mathrm{RT}}\left(\bf r,\omega\right)$
$\displaystyle=\frac{\braket{\hat{\bf{E}}^{-}\left(\bf
r,\omega\right)\hat{\bf{E}}^{+}\left(\bf{r},\omega\right)}}{\braket{\hat{\bf{E}}_{\text{p}}^{-}\left(\bf{r},\omega\right)\hat{\bf{E}}_{\text{p}}^{+}\left(\bf{r},\omega\right)}},$
(2a) $\displaystyle\text{I}_{\mathrm{RF}}\left(\bf{r},\omega\right)$
$\displaystyle=\frac{\braket{\hat{\bf{E}}^{-}_{\text{scat}}\left(\bf{r},\omega\right)\hat{\bf{E}}^{+}_{\text{scat}}\left(\bf{r},\omega\right)}}{\braket{\hat{\bf{E}}_{\text{p}}^{-}\left(\bf{r},\omega\right)\hat{\bf{E}}_{\text{p}}^{+}\left(\bf{r},\omega\right)}}.$
(2b)

or the time-resolved second-order correlation function,

$g^{(2)}(t,\tau)=\frac{\langle\hat{\bf{E}}^{-}(t)\hat{\bf{E}}^{-}(t+\tau)\hat{\bf{E}}^{+}(t+\tau)\hat{\bf{E}}^{+}(t)\rangle}{\langle\hat{\bf{E}}^{-}(t)\hat{\bf{E}}^{+}(t)\rangle^{2}}.$
(3)

In these equations, the electric-field operator can be expressed either in
frequency or time, using the standard Fourier transform relations. In either
case, the total field operator is comprised of a positive and negative
frequency component $\hat{\bf{E}}=\hat{\bf{E}}^{+}+\hat{\bf{E}}^{-}$.
Furthermore, each component of the total field can be written as a sum of the
excitation field $\hat{\bf{E}}_{\text{p}}^{\pm}$ and the field scattered (in
the low-excitation regime, as is the case in this work) or emitted (for strong
excitations) from the emitter $\hat{\bf{E}}^{\pm}_{\text{scat}}$. In the case
of RF excitation, only the scattered field arrives at the detection point.

###  II.2 Waveguide-coupled emitters

Of the different approaches that model a system of coupled quantum emitters
[25, 26, 19, 20, 27], we employ a Green’s tensor formalism developed to
describe dipolar emitters in realistic nanophotonic systems that are
dispersive, dissipative and accounting for pure dephasing of the embedded
emitters [28]. Previously, this formalism has been employed to model the
coupling of two quantum emitters via a plasmonic channel [18, 29, 24, 30], and
extended to systems of $N$ emitters coupled through a nanophotonic channel
[31, 32, 19, 20], in all cases assuming perfect coherence. We extend this
theory to include the effect of pure-dephasing for coupled emitters,
summarizing the basic principles of the theory and the way in which we account
for this decoherence.

The single emitter dipole projected Green’s function,
$\mathrm{G}(\textbf{r},\textbf{r}_{1})=\frac{i}{2}\Gamma_{1}\beta_{1}e^{ik(z-z_{1})}$
[19], can be generalized to two identical, coupled emitters in a one-
dimensional nanobeam waveguide such as is sketched in Fig. 1a and can be
written as

$\displaystyle\mathrm{G}(\bf{r}_{1},\bf{r}_{2},\omega_{\text{p}})$
$\displaystyle=i\Gamma_{0}\sqrt{\beta_{1}\beta_{2}}\frac{\hbar}{2\mu_{0}\omega_{p}^{2}}e^{ik|\Delta
z_{12}|},$ (4)

where $\omega_{p}$ is the optical angular frequency, $\mu_{0}$ is the
permeability of free-space and $k=2\pi/\lambda_{\mathrm{wg}}$ is the
wavenumber of the optical mode, whose effective wavelength is
$\lambda_{\mathrm{wg}}$. Due to the quasi-1D nature of the waveguide, the
phase of the Green’s function only depends on the emitter separation $\Delta
z$, while the emission rate $\Gamma_{0}$ and emitter-waveguide coupling
efficiencies $\beta_{i}=\Gamma_{\mathrm{wg,}i}/\Gamma_{0}$ depend implicitly
on the positions of the emitters $\bf{r}_{1}$ and $\bf{r}_{2}$.

Knowledge of this Green’s function allows us to calculate both the electric
field operator and the dispersive and dissipative coupling terms. If we assume
transverse ($\hat{\bf{y}}$-oriented) transition dipole moments $\bf{d}_{1}$
and $\bf{d}_{2}$, then the positive electric field operator isgiven by [33,
34],

$\displaystyle\hat{\bf{E}}^{+}_{\text{y}}(\bf{r})$
$\displaystyle=e^{ikz}(\hat{\bf{E}}^{+}_{\text{p0,y}}+\frac{i\hbar\Gamma_{1}\beta_{1}e^{-ikz_{1}}}{2}\hat{\sigma}_{ge}^{1}$
(5)
$\displaystyle+\frac{i\hbar\Gamma_{2}\beta_{2}e^{-ikz_{2}}}{2}\hat{\sigma}_{ge}^{2}),$

with the negative component given by the Hermitian conjugate of this
expression. Here, $\hat{\sigma}_{ge}^{i}=\ket{g_{i}}\bra{e_{i}}$ is the
lowering operator for emitter $i$. As is outlined in Appendix A, these emitter
operators are analysed in the full system Hamiltonian and contain information
about the emitter coupling rate, as well as the effects of dephasing and
excitation strength $\Omega_{p,i}$. We identify the first term in Eq. 5 as the
incident pump field (which is only present in the RT configuration), while the
following two terms represent the fields scattered from emitters 1 and 2,
respectively. In RF and resonant reflection (RR), the same scattered field
terms are used, except in RF the relative phase of the pump fields can be
tuned independently.

Within this formalism, the dispersive and dissipative coupling terms as seen
in Fig. 1 are, [19, 34, 13, 30]

$\displaystyle J_{12}$
$\displaystyle=\frac{\mu_{0}\omega_{p}^{2}}{\hbar}\bf{d}^{*}_{1}\cdot\text{Re}[\bf{G}(\bf{r}_{1},\bf{r}_{2},\omega_{p})]\cdot\bf{d}_{2}$
(6)
$\displaystyle=\frac{1}{2}\Gamma_{0}\sqrt{\beta_{1}\beta_{2}}\text{sin}(k\Delta
z_{12}),$

and,

$\displaystyle\Gamma_{12}$
$\displaystyle=\frac{2\mu_{0}\omega_{p}^{2}}{\hbar}\bf{d}^{*}_{1}\cdot\text{Im}[\bf{G}(\bf{r}_{1},\bf{r}_{2},\omega_{p})]\cdot\bf{d}_{2}$
(7) $\displaystyle=\Gamma_{0}\sqrt{\beta_{1}\beta_{2}}\text{cos}(k\Delta
z_{12}),$

respectively, where we have used
$\mathrm{G}(\bf{r}_{1},\bf{r}_{2},\omega_{\text{p}})=\bf{d}^{*}_{1}\cdot\bf{G}(\bf{r}_{1},\bf{r}_{2},\omega_{p})\cdot\bf{d}_{2}$.
As the distance between the emitters varies, the coupling changes between
being dispersive and dissipative, where the latter results in the long-lived
sub-radiant states. Note that in more complex nanophotonic systems, such as
the photonic crystal waveguides studied in Sec. II.4, the Green’s tensor is
much more complex and Eqs. (6) and (7) are modified.

Figure 2:  Ideal system response in the low-power limit
$\Omega_{p,i}=10^{-7}$. The RF (left axis) and RT (right axis) for a) a single
emitter in a waveguide, and two emitters separated by an effective phase of c)
$\pi/2$, e) $2\pi$ and g) $35\pi/18$. $g^{(2)}(\tau)$ of the corresponding RF
on resonance $\left(\delta=0\right)$ are shown in b) - h). c) and d)
correspond to ideal dispersive coupling where the energy levels are shifted
but the decay rate is unaltered, while e) and f) represent ideal dissipative
coupling where only the decay rate is altered. In the latter case, only the
super-radiant state is populated, as reflected by the rapid rise of the
correlation function. Slightly changing the relative-phase, as shown in g) and
h), enables the population of the sub-radiant state, whose decay rate is now
finite.

Using Eqs. (6) and (7), we solve the Master’s Equation of the full system
Hamiltonian and Lindblad superoperator which includes the pure-dephasing term,
as detailed in Appendix A. We recover the atomic operators
$\hat{\sigma}_{ge}^{1}$ and $\hat{\sigma}_{ge}^{2}$ as required for Eq. (5) to
determine the population of the different states. We extract both the steady-
state values, as required for the intensity calculations, and time-dependent
expectation values of atomic operator products, which is required for the
photon-correlation function $g^{(2)}(t,\tau)$. The time dependence is obtained
by the use of the quantum regression theorem [35, 36, 37].

As a baseline, we first consider the ideal case where there is no dephasing
(i.e. $\Gamma_{\mathrm{dp}}=0$), unity radiative coupling $\beta_{i}=1$,
$i=1,2$ and at low pumping strength $\Omega_{p,i}=10^{-7}$. The results in
Fig. 2 are compared to results from the literature [31, 20]. The case of RF of
a single emitter driven by a continuous-wave field is shown in Fig. 2a (left
axis) and the resonance peak corresponds to single photons being emitted into
the waveguide. Furthermore, the RT spectral response (right axis) displays a
dip due to the destructive interference between the forwardly emitted photons
and the incoming pump field [3]. The photon correlation function
$(g^{(2)}(\tau))$ displays the well-known anti-bunching behavior of a single-
photon emitter, cf. Fig. 2b, with an internal recovery-time of the emitter
lifetime ($\Gamma_{0}^{-1}$).

When the effective phase is $k\Delta z_{12}=(m+1/2)\pi$ or $k\Delta
z_{12}=p\pi$ for $m,p$ integers, the coupling is purely dispersive or
dissipative, as shown in Figs. 2c and e, respectively. For an ideal dispersive
coupling, $J_{12}=\Gamma_{0}$ and $\Gamma_{12}=0$, meaning that the energy
degeneracy of the two transitions is lifted. In RT, this results in the
appearance of a broad resonance, where the individual dips are not
individually resolved. In contrast, the RF spectrum displays two clear peaks,
and the corresponding $g^{(2)}(\tau)$ measured between them at $\delta=0$
(Fig. 2d) is oscillatory due to the quantum interference between the two
transitions. For RF, the emitters are excited with a relative phase difference
corresponding to the separation, $k\Delta z_{12}$.

For ideal dissipative coupling (Fig. 2e and f) $J_{12}=0$ and
$\Gamma_{12}=\Gamma_{0}$, meaning that both the sub-radiant $\ket{a}$ and
super-radiant $\ket{s}$ transitions are available. However, both RT and RF
spectra are determined purely by the super-radiant transitions that are
broadened compared to the single emitter spectra of a). This is confirmed by
the $g^{(2)}(\tau)$ trace, which evolves more rapidly. The perfect sub-radiant
state $(\Gamma_{0}-\Gamma_{12}=0)$ cannot emit, and therefore does not
contribute to the output spectra.

Sub-radiant behaviour may be observed when deviating from the situation of
ideal dissipative coupling ($k\Delta z_{12}=2\pi$). Figures 2g and h show
$k\Delta z_{12}=35\pi/18$, such that the sub-radiant state has a finite
lifetime. In this case, we observe a sharp, slightly detuned feature
corresponding to the sub-radiant state in both RF and RT co-existing with the
broadened super-radiant resonance. This interpretation is confirmed by the
$g^{(2)}(\tau)$ trace, where oscillations due to interference between the two
states decay on time scales of 100’s of $\Gamma_{0}^{-1}$. Further details on
the dependence of the sub-radiant peak position and width, as well as the
decay rate as seen in the $g^{(2)}(\tau)$ calculations, on the relative phase
between the emitters are given in Appendix B.

###  II.3  Role of dephasing and imperfect radiative coupling

We now consider the case where dephasing is present. Fig. 3 shows frequency-
dependent RT when $\Gamma_{\mathrm{deph}}$ ranges between 0 and
$0.5\Gamma_{0}$ (see Appendix B for the corresponding reflections).

Figure 3:  sub-radiant features in intensity RT in the presence of dephasing.
a) $I_{\mathrm{RT}}$ spectra for different $\Gamma_{\mathrm{deph}}$ for
$k\Delta z_{12}=35\pi/18$. The sharp sub-radiant feature, of magnitude $\Delta
T_{\mathrm{sub}}$, rapidly vanishes as $\Gamma_{\mathrm{deph}}$ increases, as
shown in b).  Figure 4:  The observed $g^{(2)}(\tau)$ of the coupled emitter
system when the emitter separation is $\Delta z_{12}=\pi/k$ and the emitters
are probed in RF, where the pump fields are in phase, resulting in the
excitation of the symmetric (sub-radiant) state. The observed antibunching has
a slow rise time that varies with a) dephasing $\Gamma_{deph}$ and b) $\beta$
($=\beta_{1}=\beta_{2}$). For comparison, the ideal, single-emitter
$g^{(2)}(\tau)$ trace (c.f. Fig. 2b) is shown by the dashed curve. c) 2D map
of the rise time as a function of $\Gamma_{deph}$ and $\beta$.

The super-radiant component of the RT spectra, shown in Fig. 3a, broadens and
becomes shallower as $\Gamma_{\mathrm{deph}}$ increases, as compared to the
single emitter case. Similarly, the sharp peak indicating the existence of a
sub-radiant state, $\Delta T_{\mathrm{sub}}$, rapidly vanishes as
$\Gamma_{\mathrm{deph}}$ increases, as can be seen in Fig. 3b, where we plot
the magnitude of the sub-radiant feature, $\Delta T_{\mathrm{sub}}$, as a
function of $\Gamma_{\mathrm{deph}}$. For state-of-the-art values of
$\Gamma_{\mathrm{deph}}=0.016\Gamma_{0}$ [38], we expect $\Delta
T_{\mathrm{sub}}=0.32$. Note that residual spectral diffusion, due to a slow
drift present in many experimental implementations, will reduce the sub-
radiant feature further.

In contrast, and as we observed in Fig. 2, the signature of sub-radiance in
the time dynamics, such as the photon-correlation function $g^{(2)}(\tau)$, is
pronounced and leads to a slow modulation of the trace. Furthermore, in the RF
configuration, the individual transitions can be addressed by tuning the
relative phase of the beams that excite each emitter, as discussed above.
Figure 4 shows the time-dependent $g^{(2)}(\tau)$ as a function of a)
$\Gamma_{\mathrm{deph}}$ (for almost ideal coupling, $\beta=0.99$) and b)
$\beta$ (when $\Gamma_{\mathrm{deph}}=0$), for two emitters separated by
$\Delta z_{12}=\pi/k$ and where the emitters are pumped in phase, such that
the sub-radiant state $\ket{s}$ is excited. (Note that the x-axis does not
extend all the way to zero).

Figure 5:  The observed population dynamics of the super- and sub-radiant
states, starting with both emitters initially in the ground state, with a
separation $\Delta z_{12}=\pi/k$. a) When the emitters are excited in phase
($\Omega_{p,i}=0.01$), the population of the symmetric (sub-radiant) state
reaches a steady state population of 0.5 and the super-radiant state is not
populated. b) By introducing pure dephasing, the super-radiant state is also
populated as seen in the plotted steady-state values. Inset: maximum steady-
state population (given by colorbar) of the super-radiant state for different
$\Gamma_{deph}$ and $\Omega_{p,i}$.

As the imperfections of the coupled system increase, the rise-time shortens
(see Appendix C for definition). For example, when increasing
$\Gamma_{\mathrm{deph}}$ from 0 to $0.01\Gamma_{0}$ to $0.1\Gamma_{0}$, the
rise-time changes from $100\Gamma_{0}^{-1}$ to $38.9\Gamma_{0}^{-1}$ to $6.2$
$\Gamma_{0}^{-1}$ ($\beta$=0.99). Similarly, changing $\beta$ from 0.99 to 0.9
to 0.7 changes the rise-time from $100\Gamma_{0}^{-1}$ to $10\Gamma_{0}^{-1}$
to $3\Gamma_{0}^{-1}$ ($\Gamma_{\mathrm{deph}}$=0). Encouragingly, for all of
these values, and even for $\Gamma_{deph}=0.3\Gamma_{0}$ or $\beta=0.5$, the
rise-time is significantly slower than the single-emitter response (c.f. white
curve in Fig. 4)), demonstrating that sub-radiance could be measured from the
dynamics of the coupled system. In fact, as shown in Fig. 4c, it is only in
the region where $\Gamma_{\mathrm{deph}}\approx 0.275\Gamma_{0}$ and
$\beta\approx 0.4$ that the rise-time approaches $\Gamma_{0}^{-1}$.

To further illuminate the role played by pure dephasing in the dynamics of the
coupled system, we calculate the population of symmetric and asymmetric
states, $\rho_{\mathrm{ss}}$ and $\rho_{\mathrm{aa}}$ (see Appendix A for
details), showing the results in Fig. 5. Here, the emitters are again
separated by $\Delta z_{12}=\pi/k$ and are initially in the ground state. At
time $t=0$ the pump is turned on and excite both emitters with the same local
phase. As expected, in the case of no dephasing and at weak pumping
$\Omega_{p,i}=0.01\Gamma_{0}$, we observe that the population of the sub-
radiant state oscillates to a steady-state value $\rho_{\mathrm{sub}}=0.5$
while the super-radiant state remains at $\rho_{\mathrm{sup}}=0$.

A rapid loss in steady-state population in the presence of pure dephasing is
seen in Fig. 5b, where $\rho_{\mathrm{sub}}$ drops to 0.04 at
$\Gamma_{\mathrm{deph}}=0.1\Gamma_{0}$. Interestingly, increasing the
excitation intensity can help to mitigate the effects of dephasing. That is,
the imperfections of the coupled system can, to some degree, be mitigated
through a careful selection of the excitation scheme.

A further interesting consequence of dephasing is observed in Fig. 5b; even as
$\rho_{\mathrm{sub}}$ decreases due to increased dephasing,
$\rho_{\mathrm{sup}}$ increases. That is, dephasing, introduces a pathway
between $\ket{a}$ and $\ket{s}$, allowing mixing between the two states. We
observe that, for each excitation power, $\rho_{\mathrm{sup}}$ peaks as
dephasing increases and then gradually decreases. This demonstrates the
balance between the increased mixing rate and loss of coherence introduced by
the dephasing. The optimal excitation power for a given dephasing, as required
to maximize $\rho_{\mathrm{sup}}$, is given in the inset. In the opposite
scenario (not shown), where the system is initially pumped into the super-
radiant state, we would then expect dephasing to provide a pathway to populate
the sub-radiant state.

###  II.4  Photonic crystal waveguides

While many quantum photonic experiments employ the standard nanobeam
waveguides described in the preceding sections, photonic-crystal waveguides
(PhCWs) offer higher radiative coupling efficiencies [39]. The complex
nanoscopic structure of PhCWs modifies the flow of light and creates spectral
regions where slow-light (i.e. with large group indices, $n_{g}$) modes
enhance quantum light-matter interactions [1, 40]. The electromagnetic field
is spatially modified with wavelength-scale feature sizes [41]. As a
consequence, the Green’s tensor is no longer of the form given by Eq. (4) and
is obtained from numerical solutions to the field distribution emitted by a
dipole $\bf{d}$ at position $\bf{r}^{\prime}$ according to [42],

$\bm{\text{{E}}}\left(\bf{r}\right)=\omega^{2}\mu\mu_{0}\bf{G}(\bf{r},\bf{r}^{\prime},\omega_{p})\cdot\bf{d}.$
(8)

Once the electric field distribution is known, the Green’s tensor can be
calculated using Eq. 8, instead of the simplified form given by Eq. (5).

We compute the $\hat{\bf{y}}$ component of the Green’s tensor numerically by
calculating the field emitted by a $\hat{\bf{y}}$-oriented dipole embedded in
a PhCW [41], showing the results in Fig. 6 for both a fast-
$\left(\text{n}_{\text{g}}=5\right)$ and slow-light
$\left(\text{n}_{\text{g}}=58\right)$ mode (see Appendix D for more details).

Figure 6:  Components of the Green’s tensor of a semiconductor photonic-
crystal waveguide. a) and b) show a spatial map of
$\text{Re}\>\bm{G}_{y}\left(\bm{r},\bm{r}^{\prime}\right)$, on the plane of
symmetry in the middle of a PhCW membrane, calculated for a dipole emitter 1
at position $\bm{r}^{\prime}$, for a fast-
$\left(\text{n}_{\text{g}}=5\right)$ and slow-light
$\left(\text{n}_{\text{g}}=58\right)$ mode. The axes are normalized to the
lattice constant $a$, the air holes are shown as the black circles. The second
emitter is located 10 unit cells away from emitter 1 in an area denoted by a
dashed rectangle. c) Zoom in on the dashed area with possible emitter 2
positions numerated. Note that in a) and b) emitter 1 is placed at position 4
within its unit cell. d) and e) show the position dependent dissipative and
dispersive coupling terms, $\Gamma_{12}$ and $J_{12}$, for two different group
index values of $n_{g}=5$ and $n_{g}=58$. Each square displays the position
dependent $\Gamma_{12}$/$J_{12}$ for emitter 2, where we consider two
positions for emitter 1 which are outlined by a blue box.

As expected, the Green’s tensor of the fast-light mode a) is confined to the
center of the waveguide, while that of the slow-light mode b) is more
delocalized. Note that these Green’s tensor maps are for emitter 1 located at
position 4 within the unit cell (c.f. Fig. 6c) at an anti-node of the mode.
These were repeated at other positions for emitter 1, including position 10
away from the mode maximum and near an air hole.

Once the Green’s tensor is known, $J_{12}$ and $\Gamma_{12}$ can be calculated
according to Eqs. (6) and (7), depending on the position chosen for emitter 2.
Figure 6 shows the results of such a calculation, for emitter 1 at positions 4
and 10 and emitter 2 at all positions, for both slow and fast-light modes,
where $J_{12}$ and $\Gamma_{12}$ are normalised to the natural decay rates of
emitters 1 and 2 ($\Gamma_{1}$ and $\Gamma_{2}$) within their respective unit
cells, which varies depending on the position. These maps reveal a complex
relationship between the dispersive and dissipative responses of the coupled
system in the PhCW that depends on the position and group velocity in a way
that is not present in the homogeneous nanobeam waveguides.

The coupling provided by the PhCW most closely resembles that of a nanobeam
waveguide when emitter 1 is placed at the anti-node of the mode, seen as
position 4 in Fig. 6c. For different positions of emitter 2
($\mathrm{n}_{\mathrm{g}}=5$), we observe that $\Gamma_{12}$ and $J_{12}$ are
out of phase, such that nearly perfect dissipative coupling can be achieved
(left-most squares in Figs. 6d and e respectively). An exemplary case is when
both emitters are placed at position 4,
$\Gamma_{12}\approx\Gamma_{1}(=\Gamma_{2})$ and $J_{12}\approx 0$, meaning
that near-perfect sub-radiance is possible, as we saw with the nanobeam
waveguide (c.f. Figs. 4 and 5). Emitter 2 may be displaced, for example to
positions 8 or 12, with only a slight reduction of the value of $\Gamma_{12}$,
verifying the robustness of PhCW. Alternatively, emitter 2 may be moved to
positions 2, 6 or 10, in which case $\Gamma_{12}\approx 0$ and $J_{12}\approx
0.6\sqrt{\Gamma_{1}\Gamma_{2}}$, which tends towards a dispersively coupled
system. A similar coupling behavior is observed when emitter 1 is placed away
from the mode-maximum (i.e. position 10), albeit with reduced amplitudes and
at different positions of emitter 2.

Interestingly, sub-radiant states may be easier to induce using slow-light
modes ($\mathrm{n}_{\mathrm{g}}=58$). As we see in Fig. 6d and e (third square
from the left), near-ideal dissipative coupling can still be achieved when
emitter 1 is placed at position 4, while emitter 2 is in position 3
($\Gamma_{12}\approx\sqrt{\Gamma_{1}\Gamma_{2}}$). Interestingly, displacing
emitter 2 to positions 2 or 4 does not significantly alter $\Gamma_{12}$, yet
it changes $J_{12}$ to nearly $\pm 0.2\sqrt{\Gamma_{1}\Gamma_{2}}$. Hence, by
operating in the slow-light regime in a PhCW, one can lift the energy
degeneracy of the super and sub-radiant modes. A similar but inverse trend is
seen when emitter 1 is placed away from the mode maximum, near the air hole at
position 10. Near-ideal sub-radiance is possible when emitter 2 is at
positions 1, 5 or 9, yet in these scenarios we predict that $J_{12}\approx
0.2\Gamma_{0}$, again suggesting that a spectral separation between the super-
and sub-radiant states is expected. These results indicate that PhCWs unlock
the power of complex nanoscopic structures by realizing all effective emitter
separations within one unit cell whilst maintaining a high coupling efficiency
[39, 41].

##  III  Conclusions

In summary, we have theoretically studied the possibility of forming, exciting
and detecting sub-radiant states using realistic (that is, imperfect) quantum
emitters in nanophotonic waveguides. We find that dephasing rapidly
deterioates signatures of sub-radiant states in spectral measurements, but
that they persist in the dynamics observed in the photon statistics. We
further showed that, even in the presence of dephasing, sub-radiant states may
be efficiently populated if the emitters can be excited independently,
although the dephasing leads to mixing between the super and sub-radiant
states.

We have also numerically calculated the Green’s tensor for a realistic
photonic-crystal waveguide. We showed that with a careful selection of the
position of the emitters, sub-radiance can be achieved while the energy
degeneracy of the super- and sub-radiant states may be lifted. This study,
combined with recent demonstrations of deterministic quantum-emitter photonic
structure integration [4, 43, 44, 12, 45] provides a route to the controllable
creation of sub-radiant states on-chip, providing a valuable resource to
future quantum technologies.

###### Acknowledgements.

The authors thank Anders Sondberg Sorensen, Sumanta Das, Matthew Foreman, and
Daryl M Beggs for insightful discussions and Alisa Javadi for supplying the
comsol simulation. The authors acknowledge financial support from Danmarks
Grundforskningsfond (DNRF 139, Hy-Q Center for Hybrid Quantum Networks) and
the EU’s Horizon 2020 research and innovation programme (grant No. 824140,
TOCHA, H2020-FETPROACT-01-2018). NR acknowledges funding from the Canadian
Foundation for Innovation (CFI) and the Natural Sciences and Engineering
Research Council of Canada (NSERC).

##  Appendix A Hamiltonian and Lindblad operators

To obtain an expression for the field operator, we evaluate the rising and
lowering operators, $\hat{\sigma}^{i}_{eg}$ and $\hat{\sigma}^{i}_{ge}$ by
turning to the master equation that describes the evolution of the atomic
density matrix $\hat{\rho}$ [24, 30, 19]:

$\displaystyle\frac{d}{dt}\hat{\rho}=-\frac{i}{\hbar}[\mathcal{H},\hat{\rho}]+\mathcal{L},$
(9)

where $\mathcal{H}$ and $\mathcal{L}$ are the Hamiltonian and Lindblad
operator respectively. In this paper, the results were calculated using the
atomic density operator for two emitters written in the following basis:

$\displaystyle\hat{\rho}=\left(\begin{array}[]{cccc}\rho_{11}&\rho_{12}&\rho_{13}&\rho_{14}\\\
\rho_{21}&\rho_{22}&\rho_{23}&\rho_{24}\\\
\rho_{31}&\rho_{32}&\rho_{33}&\rho_{34}\\\
\rho_{41}&\rho_{42}&\rho_{43}&\rho_{44}\\\
\end{array}\right)\quad\text{and}\quad\begin{array}[]{cc}\ket{1}&=\ket{g_{1}g_{2}}\\\
\ket{2}&=\ket{e_{1}e_{2}}\\\ \ket{3}&=\ket{g_{1}e_{2}}\\\
\ket{4}&=\ket{e_{1}g_{2}},\\\ \end{array}$ (18)

where the subscripts 1 and 2 refer to emitter 1 and 2 respectively. Using this
basis and some algebra, we can write out the relevant atomic operators for
emitter 1 and 2 in matrix form:

$\displaystyle\hat{\sigma}^{1}_{ge}=|g_{1}\rangle\langle
e_{1}|=\left(\begin{array}[]{cccc}0&0&0&1\\\ 0&0&0&0\\\ 0&1&0&0\\\ 0&0&0&0\\\
\end{array}\right)$ (23)
$\displaystyle\hat{\sigma}^{1}_{eg}=|e_{1}\rangle\langle
g_{1}|=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&1&0\\\ 0&0&0&0\\\ 1&0&0&0\\\
\end{array}\right)$ (28)
$\displaystyle\hat{\sigma}^{2}_{ge}=|g_{2}\rangle\langle
e_{2}|=\left(\begin{array}[]{cccc}0&0&1&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&1&0&0\\\
\end{array}\right)$ (33)
$\displaystyle\hat{\sigma}^{2}_{eg}=|e_{2}\rangle\langle
g_{2}|=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&1\\\ 1&0&0&0\\\ 0&0&0&0\\\
\end{array}\right)$ (38)
$\displaystyle\hat{\sigma}^{1}_{ee}=|e_{1}\rangle\langle
e_{1}|=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&1&0&0\\\ 0&0&0&0\\\ 0&0&0&1\\\
\end{array}\right)$ (43)
$\displaystyle\hat{\sigma}^{2}_{ee}=|e_{2}\rangle\langle
e_{2}|=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&0\\\
\end{array}\right).$ (48)

Here we have assumed that each emitter is a two-level system with a ground
state $|g\rangle$ and excited state $|e\rangle$, and the atomic operator
describes the transition between these two states.

We can obtain an effective Hamiltonian that describes our system within the
Born-Markov approximation, which is valid when the relevant Green’s function
does not change over the emitter linewidth [34]. This holds true for
waveguides, except in photonic crystal waveguides close to the band-edge of
the propagating mode [40]. For a system of two dipoles in a waveguide, the
effective Hamiltonian can be written as [19, 18, 34]:

$\displaystyle\mathcal{H}$
$\displaystyle=-\hbar\sum_{i=1,2}\Delta_{i}\hat{\sigma}^{i}_{ee}-\hbar\sum_{i,j=1,2}J_{ij}\hat{\sigma}^{i}_{eg}\hat{\sigma}^{j}_{ge}$
(50)
$\displaystyle-\hbar\sum_{i=1,2}(\hat{\Omega}_{p,i}^{*}\hat{\sigma}^{i}_{ge}+\hat{\Omega}_{p,i}\hat{\sigma}^{i}_{eg}),$

where $\Delta_{i}=\omega_{p}-\omega_{i}$ is the frequency detuning between
emitter $i$ and the probe field and
$\hat{\Omega}_{p,i}=\bm{d}_{i}^{*}\cdot\hat{\bm{E}}_{p}^{+}(\bm{r_{i}})/\hbar$
is the guided mode Rabi frequency. The Hamiltonian also includes the
dispersive coupling term $J_{12}$, which describes the energy level shifts due
to the coherent waveguide-mediated dipole-dipole coupling and can be evaluated
from the Green’s function[19, 34, 13, 30], as detailed in the main text. We
use the following Lindblad operator which contain the decay paths of the
system:

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{i,j=1,2}\frac{\Gamma_{ij}}{2}(2\hat{\sigma}^{i}_{ge}\hat{\rho}\hat{\sigma}^{j}_{eg}-\hat{\sigma}^{i}_{eg}\hat{\sigma}^{j}_{ge}\hat{\rho}-\hat{\rho}\hat{\sigma}^{i}_{eg}\hat{\sigma}^{j}_{ge})$
(51)
$\displaystyle+\sum_{i=1,2}\frac{\Gamma_{i,gg}}{2}(2\hat{\sigma}^{i}_{gg}\hat{\rho}\hat{\sigma}^{i}_{gg}-\hat{\sigma}^{i}_{gg}\hat{\sigma}^{i}_{gg}\hat{\rho}-\hat{\rho}\hat{\sigma}^{i}_{gg}\hat{\sigma}^{i}_{gg})$
$\displaystyle+\sum_{i=1,2}\frac{\Gamma_{i,ee}}{2}(2\hat{\sigma}^{i}_{ee}\hat{\rho}\hat{\sigma}^{i}_{ee}-\hat{\sigma}^{i}_{ee}\hat{\sigma}^{i}_{ee}\hat{\rho}-\hat{\rho}\hat{\sigma}^{i}_{ee}\hat{\sigma}^{i}_{ee}).$

Here $\Gamma_{ii}$ is the natural decay rate of each emitter and
$\Gamma_{i,gg}$, $\Gamma_{i,ee}$ are the pure dephasing rates of the ground
and excited states respectively. In realistic solid-state systems, coupling to
a phonon bath gives rise to pure dephasing processes which decrease the
coherence of the emitter [1].

By evaluating the master equation, we obtain a set of rate equations that
describe the evolution of the density matrix. The equations of motion of the
density operator can be solved for the steady-state case, where we assume the
system is excited by a continuous wave laser. However, in order to find the
steady state solution, i.e. $\dot{\hat{\rho}}=0$, we first re-write the
Hamiltonian and Lindblad operator in the so-called Dicke state basis [22]:

$\displaystyle\mathcal{H^{\prime}}$
$\displaystyle=\mathbf{B}^{-1}\mathcal{H}\mathbf{B}$ $\displaystyle
L^{\prime}$ $\displaystyle=\mathbf{B}^{-1}L\mathbf{B}$
$\displaystyle\hat{\rho}^{\prime}$
$\displaystyle=\mathbf{B}^{-1}\hat{\rho}\mathbf{B}$
$\displaystyle\rightarrow\dot{\hat{\rho}}^{\prime}$
$\displaystyle=-\frac{i}{\hbar}[\mathcal{H^{\prime}},\hat{\rho}^{\prime}]+\mathcal{L^{\prime}}$
(52)

Figure 7:  Sub-radiant features in intensity measurements in the presence of
dephasing. a) Reflection spectra, $I_{\mathrm{RR}}$, for different
$\Gamma_{\mathrm{deph}}$ for $k\Delta z_{12}=35\pi/18$. The sharp sub-radiant
feature, whose magnitude we label $\Delta T_{\mathrm{sub}}$, rapidly vanishes
as $\Gamma_{\mathrm{deph}}$ increases, also shown in $I_{RT}$ in b) for
different separations. Inset: the peak position follows the expected energy
shift, $J_{12}$, and its width corresponds to the modified decay rate
$\Gamma_{0}-\Gamma_{12}$. c) The corresponding RT $g^{(2)}(\tau)$ for the same
separations $\Delta z_{12}$ show a large bunching peak at zero time delay.

The new basis now consists of symmetric and asymmetric intermediate states,
which are entangled:

$\displaystyle\ket{g}$ $\displaystyle=\ket{g_{1}g_{2}}$ $\displaystyle\ket{e}$
$\displaystyle=\ket{e_{1}e_{2}}$ $\displaystyle\ket{s}$
$\displaystyle=\frac{1}{\sqrt{2}}(\ket{e_{1}g_{2}}+\ket{g_{1}e_{2}})$
$\displaystyle\ket{a}$
$\displaystyle=\frac{1}{\sqrt{2}}(\ket{e_{1}g_{2}}-\ket{g_{1}e_{2}})$ (53)

The steady state solution is hence obtained by setting
$\frac{d}{dt}\hat{\rho}^{\prime}=0$ and solving the 16 coupled equations
numerically. By numerically solving Eq. 9 to obtain the density matrix $\rho$,
the time evolution of the populations of the entangled states $\ket{s}$ and
$\ket{a}$, as shown in Fig. 5, can also be acquired.

##  Appendix B Sub-radiance: dependence on dephasing and separation

The reflection spectra, $\text{I}_{\mathrm{RR}}$, can be calculated using the
same equation as $\text{I}_{\mathrm{RT}}$ and is displayed in Fig. 7a. The
sharp sub-radiant feature vanishes even for moderate values of
$\Gamma_{\mathrm{deph}}$. The emitter separation, $\Delta z_{\mathrm{12}}$,
determines the ratio of dispersive coupling $J_{12}$ and dissipative coupling
$\Gamma_{12}$, as is displayed in Fig. 7b. The sharp sub-radiant feature is
primarily visible close to $k\Delta z=2\pi$, where the dissipative coupling
dominates and hence giving rise to maximum sub-radiance, as seen in the inset.
Similarly, a logarithmic decay at longer time scales is primarily a feature of
the sub-radiance, which is also dependent on the different emitter
separations.

##  Appendix C Rise-time fit

For Fig. 4, we fit the individual $g^{(2)}(\tau)$ of the coupled emitter
system to obtain the rise time. In order to compare it to the single-emitter
antibunching, we use the following function [46]:

$f(\tau)=1-\frac{\Gamma
e^{-\tau\Gamma/2}-\frac{\Gamma}{2}e^{-\tau\Gamma}}{\Gamma/2},$ (54)

where the rise time is given by the fit parameter $\Gamma$.

##  Appendix D Photonic crystal waveguide comsol simulation

We performed a full 3D numerical simulation in Comsol Multiphysics based on
the work in reference [41] to calculate the field emitted by a dipole over all
space. This allows us to use Eq. 8 to calculate the corresponding Green’s
tensor component for the PhCW.

The PhCW details that we simulate is as follows: The refractive index was
chosen to be 3.5, corresponding to the value for GaAs. The full simulation
volume spans 34 unit cells in the x-direction and 17 rows of air holes (with
refractive index 1) in the y-direction, where the unit cell and hole sizes are
240 nm and 160 nm respectively. The simulation method is briefly described as
follows: First an eigenvalue calculation is performed to determine the
eigenfrequency, eigenvector of the primary guided mode and the group index
$\mathrm{n}_{\mathrm{g}}$, from which the correct boundaries are obtained. A
finite element frequency domain simulation is then carried out with a dipole
in the PhCW using the correct boundary conditions. The simulation is repeated
for $\mathrm{n}_{\mathrm{g}}=5$ and $58$ and for different dipole positions
within the unit cell. We consider the dispersion relation:

$\mathrm{n}_{\mathrm{g}}=c/\mathrm{v}_{\mathrm{g}}=c/(\delta\omega/\delta k),$
(55)

where $\omega$ is the normalized frequency, k is the wave vector, $c$ is the
speed of light in vacuum and $\mathrm{v}_{\mathrm{g}}$ is the group velocity
that can be obtained by the slope of the guide mode. Finally, the dispersive
and dissipative coupling terms are normalised to the individual decay rates of
emitters 1 and 2 ($\Gamma_{1}$ and $\Gamma_{2}$) positioned in their
respective unit cell.

## References

  * Lodahl  _et al._ [2015] P. Lodahl, S. Mahmoodian, and S. Stobbe, Interfacing single photons and single quantum dots with photonic nanostructures, Rev. Mod. Phys. 87, 347 (2015).
  * Schmidgall  _et al._ [2018] E. R. Schmidgall, S. Chakravarthi, M. Gould, I. R. Christen, K. Hestroffer, F. Hatami, and K.-M. C. Fu, Frequency control of single quantum emitters in integrated photonic circuits, Nano Lett. 18, 1175 (2018).
  * Türschmann  _et al._ [2019] P. Türschmann, H. L. Jeannic, S. F. Simonsen, H. R. Haakh, S. Götzinger, V. Sandoghdar, P. Lodahl, and N. Rotenberg, Coherent nonlinear optics of quantum emitters in nanophotonic waveguides, Nanophotonics 8, 1641 (2019).
  * Sipahigil  _et al._ [2016] A. Sipahigil, R. E. Evans, D. D. Sukachev, M. J. Burek, J. Borregaard, M. K. Bhaskar, C. T. Nguyen, J. L. Pacheco, H. A. Atikian, C. Meuwly, R. M. Camacho, F. Jelezko, E. Bielejec, H. Park, M. Lončar, and M. D. Lukin, An integrated diamond nanophotonics platform for quantum optical networks, Science 354, 847 (2016).
  * Kim  _et al._ [2018] J.-H. Kim, S. Aghaeimeibodi, C. J. K. Richardson, R. P. Leavitt, and E. Waks, Super-radiant emission from quantum dots in a nanophotonic waveguide, Nano Lett. 18, 4734 (2018).
  * Grim  _et al._ [2019] J. Q. Grim, A. S. Bracker, M. Zalalutdinov, S. G. Carter, A. C. Kozen, M. Kim, C. S. Kim, J. T. Mlack, M. Yakes, B. Lee, and D. Gammon, Scalable in operando strain tuning in nanophotonic waveguides enabling three- quantum-dot superradiance, Nature Mat. 18, 963 (2019).
  * Trebbia  _et al._ [2022] J.-B. Trebbia, Q. Deplano, P. Tamarat, and B. Lounis, Tailoring the superradiant and subradiant nature of two coherently coupled quantum emitters, Nature Commun. 13, 2962 (2022).
  * Lund-Hansen  _et al._ [2008] T. Lund-Hansen, S. Stobbe, B. Julsgaard, H. Thyrrestrup, T. Sünner, M. Kamp, A. Forchel, and P. Lodahl, Experimental realization of highly efficient broadband coupling of single quantum dots to a photonic crystal waveguide, Phys. Rev. Lett. 101, 113903 (2008).
  * Türschmann  _et al._ [2017] P. Türschmann, N. Rotenberg, J. Renger, I. Harder, O. Lohse, T. Utikal, S. Götzinger, and V. Sandoghdar, Chip-based all-optical control of single molecules coherently coupled to a nanoguide, Nano Lett. 17, 4941 (2017).
  * Pennetta  _et al._ [2022] R. Pennetta, M. Blaha, A. Johnson, D. Lechner, P. Schneeweiss, J. Volz, and A. Rauschenbeutel, Collective radiative dynamics of an ensemble of cold atoms coupled to an optical waveguide, Phys. Rev. Lett. 128, 073601 (2022).
  * Burek  _et al._ [2017] M. J. Burek, C. Meuwly, R. E. Evans, M. K. Bhaskar, A. Sipahigil, S. Meesala, B. Machielse, D. D. Sukachev, C. T. Nguyen, J. L. Pacheco, E. Bielejec, M. D. Lukin, and M. Lončar, Fiber-coupled diamond quantum nanophotonic interface, Phys. Rev. Applied 8, 024026 (2017).
  * Pregnolato  _et al._ [2020] T. Pregnolato, X.-L. Chu, T. Schröder, R. Schott, A. D. Wieck, A. Ludwig, P. Lodahl, and N. Rotenberg, Deterministic positioning of nanophotonic waveguides around single self-assembled quantum dots, APL Photonics 5, 086101 (2020).
  * Türschmann [2017] P. Türschmann,  _Coherent Coupling of Single Organic Dye Molecules to Optical Nanoguides_ , Ph.D. thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany (2017).
  * Thyrrestrup  _et al._ [2018] H. Thyrrestrup, G. Kiršanske, H. L. Jeannic, T. Pregnolato, L. Zhai, L. Raahauge, L. Midolo, N. Rotenberg, A. Javadi, R. Schott, A. D. Wieck, A. Ludwig, M. C. Löbl, I. Söllner, R. J. Warburton, and P. Lodahl, Quantum optics with near-lifetime-limited quantum-dot transitions in a nanophotonic waveguide, Nano Lett. 18, 1801 (2018).
  * Gimeno-Segovia  _et al._ [2019] M. Gimeno-Segovia, T. Rudolph, and S. E. Economou, Deterministic generation of large-scale entangled photonic cluster state from interacting solid state emitters, Phys. Rev. Lett. 123, 070501 (2019).
  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * Beige  _et al._ [2000] A. Beige, D. Braun, B. Tregenna, and P. L. Knight, Quantum computing using dissipation to remain in a decoherence-free subspace, Phys. Rev. Lett. 85, 1762 (2000).
  * Dzsotjan  _et al._ [2010] D. Dzsotjan, A. S. Sørensen, and M. Fleischhauer, Quantum emitters coupled to surface plasmons of a nanowire: A green’s function approach, Phys. Rev. B 82, 075427 (2010).
  * Asenjo-Garcia  _et al._ [2017a] A. Asenjo-Garcia, J. D. Hood, D. E. Chang, and H. J. Kimble, Atom-light interactions in quasi-one-dimensional nanostructures: A green’s-function perspective, Phys. Rev. A 95, 033818 (2017a).
  * Das  _et al._ [2018] S. Das, V. E. Elfving, F. Reiter, and A. S. Sørensen, Photon scattering from a system of multilevel quantum emitters. ii. application to emitters coupled to a one-dimensional waveguide, Phys. Rev. A 97, 043838 (2018).
  * Aharonovich  _et al._ [2016] I. Aharonovich, D. Englund, and M. Toth, Solid-state single-photon emitters, Nature Photon. 10, 631 (2016).
  * Dicke [1954] R. H. Dicke, Coherence in spontaneous radiation processes, Phys. Rev. 93, 99 (1954).
  * Almutairi  _et al._ [2011] K. Almutairi, R. Tanaś, and Z. Ficek, Generating two-photon entangled states in a driven two-atom system, Phys. Rev. A 84, 013831 (2011).
  * Martín-Cano  _et al._ [2011] D. Martín-Cano, A. González-Tudela, L. Martín-Moreno, F. J. García-Vidal, C. Tejedor, and E. Moreno, Dissipation-driven generation of two-qubit entanglement mediated by plasmonic waveguides, Phys. Rev. B 84, 235306 (2011).
  * Kien  _et al._ [2005] F. L. Kien, S. D. Gupta, K. P. Nayak, and K. Hakuta, Nanofiber-mediated radiative transfer between two distant atoms, Phys. Rev. A 72, 063815 (2005).
  * Zheng and Baranger [2013] H. Zheng and H. U. Baranger, Persistent quantum beats and long-distance entanglement from waveguide-mediated interactions, Phys. Rev. Lett. 110, 113601 (2013).
  * Regidor  _et al._ [2021] S. A. Regidor, G. Crowder, H. Carmichael, and S. Hughes, Modeling quantum light-matter interactions in waveguide qed with retardation, nonlinear interactions, and a time-delayed feedback: Matrix product states versus a space-discretized waveguide model, Phys. Rev. Research 3, 023030 (2021).
  * Dung  _et al._ [2002] H. T. Dung, L. Knöll, and D.-G. Welsch, Resonant dipole-dipole interaction in the presence of dispersing and absorbing surroundings, Phys. Rev. A 66, 063810 (2002).
  * Martín-Cano  _et al._ [2010] D. Martín-Cano, L. Martín-Moreno, F. J. García-Vidal, and E. Moreno, Resonance energy transfer and superradiance mediated by plasmonic nanowaveguides, Nano Lett. 10, 3129–3134 (2010).
  * Gonzalez-Tudela  _et al._ [2011] A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, Entanglement of two qubits mediated by one-dimensional plasmonic waveguides, Phys. Rev. Lett. 106, 020501 (2011).
  * Fang  _et al._ [2014] Y.-L. Fang, H. Zheng, and H. U. Baranger, One-dimensional waveguide coupled to multiple qubits: photon-photon correlations, EPJ Quantum Technology 1, 13 (2014).
  * Haakh  _et al._ [2016] H. R. Haakh, S. Faez, and V. Sandoghdar, Polaritonic normal-mode splitting and light localization in a one-dimensional nanoguide, Phys. Rev. A 94, 053840 (2016).
  * Chen  _et al._ [2013] X.-W. Chen, V. Sandoghdar, and M. Agio, Coherent interaction of light with a metallic structure coupled to a single quantum emitter: From superabsorption to cloaking, Phys. Rev. Lett. 110, 153605 (2013).
  * Hood [2017] J. D. Hood,  _Atom-light interactions in a photonic crystal waveguide_ , Ph.D. thesis, California institute of technology, Pasadena, California, US (2017).
  * Lax [1963] M. Lax, Formal theory of quantum fluctuations from a driven state, Phys. Rev. 129, 2342 (1963).
  * Cohen-Tannoudji  _et al._ [2004] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Atom-photon interactions, basic processes and applications_ (John Wiley and Sons Inc, New York, United States of America, 2004).
  * Hobson and Riley [2011] M. P. Hobson and K. F. Riley,  _Essential mathematical methods for the physical sciences_ (Cambridge University Press, Cambridge, United Kingdom, 2011).
  * Jeannic  _et al._ [2021] H. L. Jeannic, T. Ramos, S. F. Simonsen, T. Pregnolato, Z. Liu, R. Schott, A. D. Wieck, A. Ludwig, N. Rotenberg, J. J. García-Ripoll, and P. Lodahl, Experimental reconstruction of the few-photon nonlinear scattering matrix from a single quantum dot in a nanophotonic waveguide, Phys. Rev. Lett. 126, 023603 (2021).
  * Arcari  _et al._ [2014] M. Arcari, I. Söllner, A. Javadi, S. L. Hansen, S. Mahmoodian, J. Liu, H. Thyrrestrup, E. H. Lee, J. D. Song, S. Stobbe, and P. Lodahl, Near-unity coupling efficiency of a quantum emitter to a photonic crystal waveguide, Phys. Rev. Lett. 113, 093603 (2014).
  * Hood  _et al._ [2016] J. D. Hood, A. Gobana, A. Asenjo-Garcia, M. Lua, S.-P. Yua, D. E. Chang, and H. J. Kimble, Atom–atom interactions around the band edge of a photonic crystal waveguide, PNAS 113, 10507–10512 (2016).
  * Javadi  _et al._ [2018] A. Javadi, S. Mahmoodian, I. Söllner, and P. Lodahl, Numerical modeling of the coupling efficiency of single quantum emitters in photonic-crystal waveguides, JOSA B 35, 514 (2018).
  * Novotny and Hecht [2006] L. Novotny and B. Hecht,  _Principles of nano-optics_ (Cambridge University Press, Cambridge, United Kingdom, 2006).
  * Davanco  _et al._ [2017] M. Davanco, J. Liu, L. Sapienza, C.-Z. Zhang, J. V. D. M. Cardoso, V. Verma, R. Mirin, S. W. Nam, L. Liu, and K. Srinivasan, Heterogeneous integration for on-chip quantum photonic circuits with single quantum dot devices, Nat. Commun. 8, 889 (2017).
  * Schnauber  _et al._ [2018] P. Schnauber, J. Schall, S. Bounouar, T. Höhne, S.-I. Park, G.-H. Ryu, T. Heindel, S. Burger, J.-D. Song, S. Rodt, and S. Reitzenstein, Deterministic integration of quantum dots into on-chip multimode interference beamsplitters using in situ electron beam lithography, Nano Lett. 18, 2336– (2018).
  * Liu  _et al._ [2021] S. Liu, K. Srinivasan, and J. Liu, Nanoscale positioning approaches for integrating single solid-state quantum emitters with photonic nanostructures, Laser Photonics Rev. 15, 2100223 (2021).
  * Loudon [2010] R. Loudon,  _The Quantum Theory of Light_ (Oxford University Press, Oxford, United Kingdom, 2010).
  * Asenjo-Garcia  _et al._ [2017b] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays, Phys. Rev. X 7, 031024 (2017b).
  * Sun  _et al._ [2016] S. Sun, H. Kim, G. S. Solomon, and E. Waks, A quantum phase switch between a single solid-state spin and a photon, Nature Nanotech. 11, 539 (2016).
  * Elshaari  _et al._ [2017] A. W. Elshaari, I. E. Zadeh, A. Fognini, M. E. Reimer, D. Dalacu, P. J. Poole, V. Zwiller, and K. D. Jöns, On-chip single photon filtering and multiplexing in hybrid quantum photonic circuits, Nature Comm. 8, 379 (2017).
  * Elshaari  _et al._ [2020] A. W. Elshaari, W. Pernice, K. Srinivasan, O. Benson, and V. Zwiller, Hybrid integrated quantum photonic circuits, Nature Photon. 14, 285 (2020).
  * Midolo  _et al._ [2018] L. Midolo, A. Schliesser, and A. Fiore, Nano-opto-electro-mechanical systems, Nature Nanotech. 13, 11 (2018).
  * Pelton  _et al._ [2002] M. Pelton, C. Santori, J. Vučković, B. Zhang, G. S. Solomon, J. Plant, and Y. Yamamoto, Efficient source of single photons: A single quantum dot in a micropost microcavity, Phys. Rev. Lett. 89, 233602 (2002).
  * Claudon  _et al._ [2010] J. Claudon, J. Bleuse, N. S. Malik, M. Bazin, P. Jaffrennou, N. Gregersen, C. Sauvan, P. Lalanne, and J.-M. Gérard, A highly efficient single-photon source based on a quantum dot in a photonic nanowire, Nature Photon. 4, 174 (2010).
  * Versteegh  _et al._ [2014] M. A. M. Versteegh, M. E. Reimer, K. D. Jöns, D. Dalacu, P. J. Poole, A. Gulinatti, A. Giudice, and V. Zwiller, Observation of strongly entangled photon pairs from a nanowire quantum dot, Nature Comm. 5, 5298 (2014).
  * Dousse  _et al._ [2010] A. Dousse, J. Suffczyński, A. Beveratos, O. Krebs, A. Lemaître, I. Sagnes, J. Bloch, P. Voisin, and P. Senellart, Ultrabright source of entangled photon pairs, Nature 466, 216 (2010).
  * Kuhlmann  _et al._ [2013] A. V. Kuhlmann, J. Houel, A. Ludwig, L. Greuter, D. Reuter, A. D. Wieck, M. Poggio, and R. J. Warburton, Charge noise and spin noise in a semiconductor quantum device, Nature Phys. 9, 570 (2013).
  * Borregaard  _et al._ [2019] J. Borregaard, A. S. Sørensen, and P. Lodahl, Quantum networks with deterministic spin-photon interfaces, Adv. Quantum Technol. 2, 1800091 (2019).
  * Petroff  _et al._ [2001] P. M. Petroff, A. Lorke, and A. Imamoglu, Epitaxially self-assembled quantum dots, Phys. Today 54, 46 (2001).
  * Mandel  _et al._ [2003] O. Mandel, M. Greiner, A. Widera, T. Rom, T. W. Hänsch, and I. Bloch, Controlled collisions for multi-particle entanglement of optically trapped atoms, Nature 425, 937 (2003).
  * Nägerl  _et al._ [1999] H. C. Nägerl, D. Leibfried, H. Rohde, G. Thalhammer, J. Eschner, F. Schmidt-Kaler, and R. Blatt, Laser addressing of individual ions in a linear ion trap, Phys. Rev. A 60, 145 (1999).
  * Noh and Angelakis [2017] C. Noh and D. G. Angelakis, Quantum simulations and many-body physics with light, Rep. Prog. Phys. 80, 016401 (2017).
  * Mahmoodian  _et al._ [2020] S. Mahmoodian, G. Calaj’o, D. E. Chang, K. Hammerer, and A. S. Sørensen, Dynamics of many-body photon bound states in chiral waveguide qed, Phys. Rev. X 10, 031011 (2020).
  * Paulisch  _et al._ [2016] V. Paulisch, H. J. Kimble, and A. González-Tudela, Universal quantum computation in waveguide qed using decoherence free subspaces, New J. Phys. 18, 043041 (2016).
  * Hettich  _et al._ [2002] C. Hettich, C. Schmitt, J. Zitzmann, S. Kühn, I. Gerhardt, and V. Sandoghdar, Nanometer resolution and coherent optical dipole coupling of two individual molecules, Science 298, 385 (2002).


# Multi-year Investment Modelling in Energy Systems

Diego A. Tejada-Arango

(July 2023)

##  1 Nomenclature

###  1.1 Indices

Name |  Description
---|---
$y,m,\mu$ |  years
$k$ |  representative periods
$t$ |  time steps

###  1.2 Sets

Name |  Description
---|---
$\mathcal{M}$ |  Milestone years

###  1.3 Parameters

Name |  Description | Units
---|---|---
$R$ |  Social Discount Rate | [%]
$WACC_{y}$ |  Weighted Average Cost of Capital at year $y$ | [%]
$LT$ |  Lifetime of the asset | [years]
$Y$ |  Total number of years in the time horizon | [years]
$C_{y}^{T}$ |  Total investment cost at year $y$ (i.e., overnight cost) | [€/MW]
$C_{y}^{A}$ |  Annualized investment cost at year $y$ | [€/MW/year]
$C^{I}$ |  Total investment cost | [€]
$C^{O}$ |  Total operational cost | [€]
$SV_{y}$ |  Salvage value | [€]
$SVP_{y}$ |  Salvage value percentage | [%]
$C_{y}^{op}$ |  Operational cost | [€/MWh]
$W_{yk}^{op}$ |  Operational cost weight of representative $k$ in year $y$ | [h]
$W_{y}^{I}$ |  Milestone year weight | [years]
$W_{my}^{O}$ |  Weight of milestone year $m$ on year $y$ | [0-1]
$W_{m\mu}^{M}$ |  Weight of milestone year $m$ on operational year $\mu$ | [0-1]
$W_{my\mu}^{Y}$ |  Weight of milestone year $m$ on year $y$ and operational year $\mu$ | [0-1]

###  1.4 Variables

Name |  Description | Units
---|---|---
$x_{y}$ |  capacity investment at year $y$ | [MW]
$p_{ykt}$ |  production at period $t$, representative $k$, and year $y$ | [MW]
$p_{mykt}$ |  production at period $t$, representative $k$, operational year $y$, from the investment made in year $m$ | [MW]

##  2 Introduction

This document summarises the main multi-year investment modelling approaches
in energy planning models. Therefore, here we will go from a simple (basic)
formulation to a more complex (general) one to understand different levels of
detail, including examples to make more accessible the understanding of the
concepts.

In the literature, there are two approaches to managing multi-horizon
investments:

  * •

Total Cost: Using the total cost (a.k.a. overnight cost) for the investment
year.

$C_{y}^{T}=\sum_{j=y}^{y+LT-1}\frac{1}{(1+WACC_{y})^{j-y}}\cdot
C_{j}^{A}\quad\forall y$ (1)

  * •

Annualized Cost: Distributing the total cost over the years the asset is
active, which includes all the years following the build year until its
lifetime.

$C_{y}^{A}=\frac{WACC_{y}}{(1+WACC_{y})\cdot\left(1-\frac{1}{(1+WACC_{y})^{LT}}\right)}\cdot
C_{y}^{T}\quad\forall y$ (2)

It’s important to note that equation (1) has infinite solutions for defining
the annuities ($C_{j}^{A}$). One standard method is to assume that the
annuities remain constant throughout the lifetime, as illustrated in equation
(2). In addition, equation (2) does not discount the first year of the
investment, which is explained in [1]. Another common approach is discounting
the first year and using equation (3) instead. Depending on the assumptions
made, either equation (2) or (3) can be used to calculate the annualized cost
from the total cost.

$C_{y}^{A}=\frac{WACC_{y}}{1-\frac{1}{(1+WACC_{y})^{LT}}}\cdot
C_{y}^{T}\quad\forall y$ (3)

Let us consider the following data as an example for the first year ($y=0$):

  * •

Total capital cost $C_{0}^{T}=100$ [€/MW]

  * •

The $WACC_{0}=2$ %

  * •

The economic lifetime $LT=5$ years

Therefore, the annualized cost for each case is:

  * •

Annualized capital cost using equation (2): $C_{0}^{I}=20.80$ [€/MW]

  * •

Annualized capital cost using equation (3): $C_{0}^{I}=25.05$ [€/MW]

Although the values seem different, they will represent the total investment
cost if correctly discounted for each case. For instance, using the annualized
value from equation (2):

$C_{0}^{T}=100=20.80\cdot\left(1+\frac{1}{(1+0.02)^{1}}+\frac{1}{(1+0.02)^{2}}+\frac{1}{(1+0.02)^{3}}+\frac{1}{(1+0.02)^{4}}\right)$

Whereas using the annualized value from equation (3) the calculation is:

$C_{0}^{T}=100=25.05\cdot\left(\frac{1}{(1+0.02)^{1}}+\frac{1}{(1+0.02)^{2}}+\frac{1}{(1+0.02)^{3}}+\frac{1}{(1+0.02)^{4}}+\frac{1}{(1+0.02)^{5}}\right)$

Therefore, using equation (3) means that equation (1) changes slightly to
represent the total cost adequately. For the remainder of this document, we
will utilize equations (1) and (2) to continue our discussion.

The reference [2] provides a concise summary of the approach utilized in
popular planning models within the energy sector, including TIMES, DIMENSION,
SWITCH, and PRIMES. In addition, it suggests that the methods are only
comparable under the condition that $Y\rightarrow\infty$. Otherwise, this
could lead to different interpretations of the costs in the objective
function. However, this report demonstrates that with careful execution, the
objective function can be made identical in both techniques.

##  3 Basic formulations

Here we will go from a simple (basic) formulation to a more complex (detailed)
one for each method. Equation (4) shows a general formulation of the objective
function for a multi-year investment planning model, while equation (5) shows
the relationship between the variables where the investment decisions limit
generation production considering its lifetime.

$\min_{x_{y},p_{ykt}}\quad C^{I}+C^{O}$ (4)

$\displaystyle s.t.$ (5) $\displaystyle
p_{ykt}\leq\sum_{j=\max(y-LT+1,0)}^{y}x_{j}\quad\forall y\forall k\forall t$

Where $C^{I}$ refers to the investment cost and $C^{O}$ to the operational
cost. A general definition for $C^{O}$ is shown in equation (6):

$C^{O}=\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\cdot
C_{y}^{op}\sum_{k}W_{yk}^{op}\cdot\sum_{t}p_{ykt}\\\ $ (6)

As for the investment cost component in the objective function $C^{I}$, this
document will focus on obtaining this component in both approaches: the total
and annualized investment costs. There are additional costs that have not been
included in this objective function for the sake of simplicity, such as fixed
costs, decommissioning costs, and costs for mothballing.

###  3.1 Total Investment Cost

Equation (7) shows the total investment cost through the modelled years,
considering the social discount rate $R$.

$\displaystyle C^{I}=\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\cdot C_{y}^{T}x_{y}$
(7) $\displaystyle C^{I}=C_{0}^{T}x_{0}+\frac{1}{(1+R)^{1}}\cdot
C_{1}^{T}x_{1}+\frac{1}{(1+R)^{2}}\cdot C_{2}^{T}x_{2}+\cdots$

###  3.2 Annualized Investment Cost

Equation (8) is obtained by replacing (1) in (7). Bear in mind that the value
of the parameter $C_{y}^{A}$ is determined by (2).

$\displaystyle C^{I}=$ $\displaystyle\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\cdot
x_{y}\cdot\sum_{j=y}^{\min(y+LT-1,Y-1)}\frac{1}{(1+WACC_{y})^{j-y}}\cdot
C_{y}^{A}$ (8) $\displaystyle C^{I}=$ $\displaystyle
C_{0}^{A}\cdot\left(1+\frac{1}{(1+WACC_{0})^{1}}+\frac{1}{(1+WACC_{0})^{2}}+\cdots\right)\cdot
x_{0}$ $\displaystyle+\frac{1}{(1+R)^{1}}\cdot
C_{1}^{A}\cdot\left(1+\frac{1}{(1+WACC_{1})^{1}}+\frac{1}{(1+WACC_{1})^{2}}+\cdots\right)\cdot
x_{1}$ $\displaystyle+\frac{1}{(1+R)^{2}}\cdot
C_{2}^{A}\cdot\left(1+\frac{1}{(1+WACC_{2})^{1}}+\frac{1}{(1+WACC_{2})^{2}}+\cdots\right)\cdot
x_{2}+\cdots$

##  4 Analysis

As mentioned in Section 2, equations (7) and (8) are equivalent when
considering yearly investment and $Y\rightarrow\infty$. This statement can be
checked by considering the following relationships:

$\displaystyle C_{0}^{T}=$ $\displaystyle
C_{0}^{A}\cdot\left(1+\frac{1}{(1+WACC_{0})^{1}}+\frac{1}{(1+WACC_{0})^{2}}+\cdots\right)$
$\displaystyle C_{1}^{T}=$ $\displaystyle
C_{1}^{A}\cdot\left(1+\frac{1}{(1+WACC_{1})^{1}}+\frac{1}{(1+WACC_{1})^{2}}+\cdots\right)$
$\displaystyle C_{2}^{T}=$ $\displaystyle
C_{2}^{A}\cdot\left(1+\frac{1}{(1+WACC_{2})^{1}}+\frac{1}{(1+WACC_{2})^{2}}+\cdots\right)$
$\displaystyle\vdots$ $\displaystyle\vdots$

It is also essential to highlight that no matter the values of the social
discount rate ($R$) or the technology-specific discount rate ($WACC$), the
total investment cost ($C^{I}$) in the objective function is the same. This
affirmation holds if you obtain the annualized value using equation (2) and
the same $WACC$ value in the equation (8). If this isn’t the situation,
there’s a possibility of having different total values in the objective
function. In that situation, unifying the $WACC$ and using the same value for
both equations is better to avoid errors and confusion.

Finally, the $Y\rightarrow\infty$ assumption is not realistic to apply. The
usual use case considers that modelled years do not cover the entire
investment lifetime and use milestone years (or representative years) to
reduce the computational burden. These two cases are discussed in the
following sections.

###  4.1 Having Modelled Years $<$ Lifetime

First, Let’s analyze the case of years that do not cover the entire investment
lifetime but still use yearly investment decisions. For instance, if the last
modelled year ($Y$) is 2050 and we have an investment made in 2040 with a
lifetime of 15 years, then the last four years will not be considered in the
optimisation. In this situation:

  * •

The annualized cost approach in (8) will only recover the annualized
investment costs from 2040 to 2050, and therefore it assumes that the non-
modelled years (i.e., 2051 to 2054) will be recovered independently from the
modelled years and are not part of the optimization. Under this assumption,
there is no need for changing the equation (8) due to modelling years below
the asset’s lifetime.

  * •

The total investment cost approach has a problem using the equation (7) in the
objective function. It will consider the total cost of the investment but will
miss the operational costs from the non-modelled years to compensate for the
total investment cost. In our example, the operational cost from 2040 to 2050
must compensate for the total investment cost of a 15-year lifetime asset. To
avoid this situation, the standard approach is to consider the salvage value
($SV$).

####  4.1.1 The Salvage value

The salvage value is what a company expects to get by selling or disassembling
an asset at the end of its useful life, which will be the last modelled year
in our case for these types of models. As in the annualized cost approach, one
common assumption is that the salvage value is equivalent to the cost outside
the modelled years (i.e., the non-modelled years will be recovered
independently of the optimization problem). Therefore, the salvage value that
satisfies this assumption can be determined from the annualized cost value,
see equation (9).

$SV_{y}=C_{y}^{A}\cdot\sum_{j=Y+1}^{y+LT-1}\frac{1}{(1+WACC_{y})^{j-y}}\quad\forall
y\in[0,Y]$ (9)

The salvage value $SV$ will depend on the year $y$, the last modelled year
$Y$, the annualized value $C_{y}^{A}$, the lifetime $LT$, and the $WACC$. For
instance, consider the following data for the first year ($y=0$):

  * •

Total capital cost $C_{0}^{T}=100$ [€/MW]

  * •

The $WACC_{0}=5$ %

  * •

The last modelled year $Y=4$

  * •

The economic lifetime $LT=8$ years

  * •

Annualized capital cost using equation (2): $C_{0}^{I}=14.74$ [€/MW]

Therefore, using equation (9), the $SV_{0}=33.01$ [€/MW] corresponds to the
remaining cost in non-modelled years: 5, 6, and 7. The $SV$ is referenced to
the first year $y=0$ and can be considered as a benefit of the investment in
the objective function, as in equations (10) and (11).

$C^{I}=\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\cdot(C_{y}^{T}-SV_{y})\cdot x_{y}$
(10)

$C^{I}=\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\cdot C_{y}^{T}(1-SVP_{y})\cdot
x_{y}$ (11)

Where,

$SVP_{y}=\frac{SV_{y}}{C_{y}^{T}}\quad\forall y$

Equations (10) and (11) are equivalent and can be used instead of equation (7)
to account for the salvage value in the total investment cost for the
objective function.

If the $SV$ is obtained using the same assumptions as the ones for the
annualized cost method, then considering the $SV$ will make both methods
equivalent. In the final year of modelling, $Y$, the $SV$ will make up for
LT-1 years of the asset. This implies that the total investment method will
only consider the cost of one year in the last modelled year, resulting in the
same value as the annualized cost approach. As a result, both methods
represent the same cost in the objective function.

###  4.2 Using Milestone Years

When studying pathways to the year 2050, it is impractical to model yearly
investment decisions on a large scale. Therefore, a common approach is to use
the so-called milestone years, where the major investment decisions will be
made (e.g., years 2030, 2035, 2040, and 2050). Each milestone year represents
part of the non-modelled years.

Consider the information in Table 1 as an example; therefore, years 1, 3, and
5 are not modelled since the milestone years for investments are years 0, 2,
and 5.

Years $y$ | 0 | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
Milestone year $m$ | 0 | - | 2 | - | - | 5
Milestone year weight $W_{m}^{I}$ | 2 | - | 3 | - | - | 1
Table 1: Milestone years

The total investment cost approach needs a weight coefficient in the
operational cost $W_{y}^{I}$ to properly account for the non-modelled years.
In this example, year 0 has a weight of 2, year 2 has a weight of 3, whereas
year 4 equals 1. Moreover, since the total investment cost is considered in
the milestone year, the total investment cost $C^{I}$ is not needed to
multiply by the coefficient $W_{y}^{I}$. Equation (12) shows the formulation,
including all these considerations:

$\displaystyle\min_{x_{m},p_{mkt}}\quad C^{I}+C^{O}$ (12) $\displaystyle
C^{I}=\sum_{m\in\mathcal{M}}\frac{1}{(1+R)^{m}}\cdot(C_{m}^{T}-SV_{m})\cdot
x_{m}$ $\displaystyle
C^{O}=\sum_{m\in\mathcal{M}}\frac{1}{(1+R)^{m}}\cdot{\color[rgb]{0,0,1}W_{m}^{I}}\cdot
C_{m}^{op}\sum_{k}W_{mk}^{op}\cdot\sum_{t}p_{mkt}$

The annualized investment cost approach does not need a weight coefficient a
priory in its formulation if one assumes that the not modelled years recover
their own cost (i.e., the same assumption as for the years after the last
modelled year). However, it might be desirable to include the weight
coefficient $W_{m}^{I}$ for this method if the milestone years in unevenly
spaced and one wants to consider that a particular milestone year has more
weight than the others. For instance, the asset has more production costs at
the beginning and less at the end of the horizon, such as the production costs
at the beginning compensate enough for investment cost. Therefore, using the
weight coefficient $W_{m}^{I}$ in this method is also usual practice; see
equation (13).

$\displaystyle\min_{x_{m},p_{mkt}}\quad C^{I}+C^{O}$ (13) $\displaystyle
C^{I}=\sum_{m\in\mathcal{M}}\frac{1}{(1+R)^{m}}\cdot
x_{m}\cdot\sum_{j=m|j\in\mathcal{M}}^{\min(m+LT-1,Y-1)}\frac{1}{(1+WACC_{m})^{j-m}}\cdot{\color[rgb]{0,0,1}W_{j}^{I}}\cdot
C_{m}^{A}$ $\displaystyle
C^{O}=\sum_{m\in\mathcal{M}}\frac{1}{(1+R)^{m}}\cdot{\color[rgb]{0,0,1}W_{m}^{I}}\cdot
C_{m}^{op}\sum_{k}W_{mk}^{op}\cdot\sum_{t}p_{mkt}$

Please take note that the calculation for investment cost in equation (13)
only considers milestone years in the inner summation, owing to
$j\in\mathcal{M}$.

Let us continue with the example in Table 1 (i.e., $\mathcal{M}=\\{0,2,5\\}$
and $W_{m}^{I}=[2,3,1]$) and assume that $LT=6$. Therefore, the investment
cost in equations in the previous equations becomes:

Equation (12) - total investment cost method:

$C^{I}={\color[rgb]{1,.5,0}C_{0}^{T}}x_{0}+\frac{1}{(1+R)^{2}}(C_{2}^{T}-SV_{2})\cdot
x_{2}+\frac{1}{(1+R)^{5}}(C_{5}^{T}-SV_{5})\cdot x_{5}$

Equation (13) - annualized investment cost method:

$\displaystyle C^{I}=$
$\displaystyle{\color[rgb]{1,.5,0}C_{0}^{A}\left(2+\frac{3}{(1+WACC_{0})^{2}}+\frac{1}{(1+WACC_{0})^{5}}\right)}\cdot
x_{0}$ $\displaystyle+\frac{1}{(1+R)^{2}}\cdot
C_{2}^{A}\left(3+\frac{1}{(1+WACC_{2})^{3}}\right)\cdot
x_{2}+\frac{1}{(1+R)^{5}}C_{5}^{A\cdot}x_{5}$

By simply inspecting the methods used, it is clear that the costs represented
in the objective function are not the same as before when yearly investment
decisions were made. For example, if we compare the first term in both costs
(highlighted in orange), we can see that the values differ according to
equations (1) and (2).

The operational cost is the same in both methods since the assumption here is
that the milestone years represent the non-modelled years using the weight
$W_{m}^{I}$. Applying the values of the example in equations (12) and (13),
the operational cost is:

$\displaystyle C^{O}=$
$\displaystyle\sum_{m\in\mathcal{M}}\frac{1}{(1+R)^{m}}\cdot{\color[rgb]{0,0,1}W_{m}^{I}}\cdot
C_{m}^{op}\sum_{k}W_{mk}^{op}\cdot\sum_{t}p_{mkt}$ $\displaystyle+2\cdot
C_{0}^{op}\sum_{k}W_{0,k}^{op}\cdot\sum_{t}p_{0,kt}$
$\displaystyle+\frac{3}{(1+R)^{2}}\cdot
C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{2,kt}$
$\displaystyle+\frac{1}{(1+R)^{5}}\cdot
C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}p_{5,kt}$

Therefore, this example assumes that milestone year 2 represents non-modelled
years 3 and 4. The issue here is that years 3 and 4 should have a lower value
of money over time than year 2, meaning that this approach overestimates the
operational cost. Furthermore, this approach presents another issue. If the
asset’s lifetime ends between two milestone years, the $C^{I}$ and $C^{O}$ in
the objective function cannot account for this and will assume the asset is
available throughout the entire range between the milestone years. For
instance, if $LT=4$ in Table (1), the asset should be available until year 3.
However, equations (2) will include them in the objective function until the
final year (i.e., $y=5$), resulting in an additional two years in the
objective function (i.e., overestimation of the operational costs again).

##  5 Proposals to Improve the Milestone Years Method

As discussed in the previous section, multiplying the milestone years by the
weight $W_{m}^{I}$ will lead to two main issues in the objective function of
both methods:

  * •

The value of money over time is not correctly represented for the non-modelled
years.

  * •

The lifetime of the assets is overestimated if the last year is between two
milestone years.

In the following sections, we propose some formulations that attempt to solve
these shortcomings in both approaches.

###  5.1 Proposal 1 - Mapping investment milestone years to all years
($W_{my}^{O}$)

The first proposal is based on adding an extra index $y$ in the objective
function that includes all the years in the time horizon (i.e., milestone and
non-modelled years). Then, a mapping between the non-modelled years and the
milestone year is needed to properly account for them. The weight of milestone
year $m$ on year $y$ ($W_{my}^{O}$) ensures this mapping among years. For the
total investment cost, we only need to change the operational cost in the
objective function as shown in equation (14). The very change applies to the
annualized investment cost approach, which will account additionally for all
the years during the asset’s lifetime in the $C^{I}$; see equation (15).

Total investment cost method:

$\displaystyle\min_{x_{m},p_{mkt}}\quad C^{I}+C^{O}$ (14) $\displaystyle
C^{I}=\sum_{m\in\mathcal{M}}\frac{1}{(1+R)^{m}}\cdot(C_{m}^{T}-SV_{m})\cdot
x_{m}$ $\displaystyle
C^{O}=\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\sum_{m\in\mathcal{M}}{\color[rgb]{0,0,1}W_{my}^{O}}\cdot
C_{m}^{op}\sum_{k}W_{mk}^{op}\cdot\sum_{t}p_{mkt}$

Annualized investment cost method:

$\displaystyle\min_{x_{m},p_{mkt}}\quad C^{I}+C^{O}$ (15) $\displaystyle
C^{I}=\sum_{m\in\mathcal{M}}\frac{1}{(1+R)^{m}}\cdot
x_{m}\cdot\sum_{y=m}^{\min(m+LT-1,Y-1)}\frac{1}{(1+WACC_{m})^{y-m}}\cdot
C_{m}^{A}$ $\displaystyle
C^{O}=\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\sum_{m\in\mathcal{M}}{\color[rgb]{0,0,1}W_{my}^{O}}\cdot
C_{m}^{op}\sum_{k}W_{mk}^{op}\cdot\sum_{t}p_{mkt}$

Let us first analyze the $C^{I}$ in both equations assuming that the $LT=5$.
Therefore, the investment cost in equations in the previous equations becomes:

Equation (14) - total investment cost method:

$C^{I}={\color[rgb]{0.5,0.5,0}C_{0}^{T}}x_{0}+\frac{1}{(1+R)^{2}}(C_{2}^{T}-SV_{2})\cdot
x_{2}+\frac{1}{(1+R)^{5}}(C_{5}^{T}-SV_{5})\cdot x_{5}$

Equation (15) - annualized investment cost method:

$\displaystyle C^{I}=$
$\displaystyle{\color[rgb]{0.5,0.5,0}C_{0}^{A}\left(1+\frac{1}{(1+WACC_{0})^{1}}+\frac{1}{(1+WACC_{0})^{2}}+\frac{1}{(1+WACC_{0})^{3}}+\frac{1}{(1+WACC_{0})^{4}}\right)}\cdot
x_{0}$
$\displaystyle+\frac{1}{(1+R)^{2}}C_{2}^{A}\left(1+\frac{1}{(1+WACC_{2})^{1}}+\frac{1}{(1+WACC_{2})^{2}}+\frac{1}{(1+WACC_{2})^{3}}\right)\cdot
x_{2}$ $\displaystyle+\frac{1}{(1+R)^{5}}C_{5}^{A}\cdot x_{5}$

By simply inspecting the methods used, it is clear that the costs represented
in the objective function are the same as before when yearly investment
decisions were made. For example, if we compare the first term in both costs
(highlighted in olive), we can see that the values are equivalent according to
equations (1) and (2). Therefore, this proposal deals with the problem of
having an asset whose lifetime ends in the middle of two milestone years for
the annualized investment cost method in the $C^{I}$.

Regarding the operational cost $C^{O}$, Table 2 shows an example of the values
in the parameter $W_{my}^{O}$. Please note that it only impacts the
operational cost component in equations (14) and (15). For instance, the
operational cost of non-modelled year 1 will be represented as half the values
of milestone years 0 and 2. Non-modelled years 3 and 4 are also a linear
combination of the adjacent milestone years 2 and 5. This example uses a
linear distribution; however, since it is an input parameter, any other type
of distribution can be used.

| Years $y$
---|---
Milestone Year $m$ |  0 |  1 |  2 |  3 |  4 |  5
0 |  1 |  1/2 |  0 |  0 |  0 |  0
2 |  0 |  1/2 |  1 |  2/3 |  1/3 |  0
5 |  0 |  0 |  0 |  1/3 |  2/3 |  1
Table 2: Weight of milestone year $m$ on year $y$ ($W_{my}^{O}$)

The parameter $W_{my}^{O}$ allows mapping of each not modelled year to a
milestone year, considering the value of money over time. Using equations (14)
and (15) and the values of the example, we get:

$\displaystyle C^{O}=$ $\displaystyle
C_{0}^{op}\sum_{k}W_{0k}^{op}\cdot\sum_{t}p_{0kt}$
$\displaystyle+\frac{1}{(1+R)^{1}}\cdot\left(\frac{1}{2}\cdot\sum_{k}W_{0k}^{op}\cdot\sum_{t}p_{0kt}+\frac{1}{2}\cdot\sum_{k}W_{2k}^{op}\cdot\sum_{t}p_{2kt}\right)$
$\displaystyle+\frac{1}{(1+R)^{2}}\cdot\left(\sum_{k}W_{2k}^{op}\cdot\sum_{t}p_{2kt}\right)$
$\displaystyle+\frac{1}{(1+R)^{3}}\cdot\left(\frac{2}{3}\cdot\sum_{k}W_{2k}^{op}\cdot\sum_{t}p_{2kt}+\frac{1}{3}\cdot\sum_{k}W_{5k}^{op}\cdot\sum_{t}p_{5kt}\right)$
$\displaystyle+\frac{1}{(1+R)^{4}}\cdot\left(\frac{1}{3}\cdot\sum_{k}W_{2k}^{op}\cdot\sum_{t}p_{2kt}+\frac{2}{3}\cdot\sum_{k}W_{5k}^{op}\cdot\sum_{t}p_{5kt}\right)$
$\displaystyle+\frac{1}{(1+R)^{5}}\cdot\left(\sum_{k}W_{5k}^{op}\cdot\sum_{t}p_{5kt}\right)$

Here we can see the main drawback of this proposal, since equations (14) and
(15) assume that the operational cost is a linear combination of the milestone
years (i.e., $m=$ 0, 2, and 5) for non-modelled years (i.e., $y=$ 1, 3, and
4), then the non-modelled years like years 3 and 4 will consider part of the
production costs from the milestone year 5 (which is a future year). There is
an improvement compared to the standard approach of having weights in the
milestone years shown in Section (4.2); however, it still might have an
overestimation of the operational costs in the objective function. The
following proposal aims to overcome this issue.

###  5.2 Proposal 2 - Mapping investment milestone years to all years and
operational milestone years ($W_{my\mu}^{Y}$ and $W_{m\mu}^{M}$)

Equations (14) and (15) in proposal 1 solve the issue with the investment cost
($C^{I}$) since the representation of the investment cost is the same in both
methods. In addition, it improved the value of money over time representation
for the operational cost. Nevertheless, there is still a concern about
accounting operational costs from future investment decisions. One solution to
that situation is shown in equation (16).

$C^{O}=\sum_{m\in\mathcal{M}}\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\sum_{\mu\in\mathcal{M}}{\color[rgb]{0,0,1}W_{my\mu}^{Y}}\cdot
C_{m}^{op}\sum_{k}W_{mk}^{op}\cdot\sum_{t}{\color[rgb]{.75,0,.25}W_{m\mu}^{M}}\cdot
p_{mkt}$ (16)

Where,

  * •

The parameter $W_{m\mu}^{M}$ represents the proportion of production assigned
to the investment made in milestone year $m$ on the operational milestone year
$\mu$.

  * •

The parameter $W_{my\mu}^{Y}$ represents the weight of investment milestone
year $m$ on year $y$ considering operational milestone year $\mu$. In other
words, how much of the operational production of year $\mu$ does one want to
consider for the year $y$ of the investment made in year $m$.

To understand the parameter $W_{m\mu}^{M}$, let us first consider in our
example that equation (5) becomes equation (17) when taking into account the
milestone years $m=$ 0, 2, and 5.

$\displaystyle p_{0kt}\leq x_{0}$ $\displaystyle\quad\forall k\forall t$ (17)
$\displaystyle p_{2kt}\leq x_{0}+$ $\displaystyle x_{2}$
$\displaystyle\quad\forall k\forall t$ $\displaystyle p_{5kt}\leq$
$\displaystyle x_{2}+x_{5}$ $\displaystyle\quad\forall k\forall t$

Equation (17) shows that operational production in year 0 is due entirely to
the investment made in year 0, while the operational production in year 2
depends on the investments made in both year 0 and year 2, and the operational
production in year 5 depends on the investments made in year 2 and year 5, but
not year 0, since the lifetime ($LT$=5). Multiple assumptions can be made here
to determine the parameter $W_{m\mu}^{M}$. A naive approach is to distribute
equally the number of investments that limits each operational year. For
instance, the operational production of year 2 can be split into halves
between the investments made in years 0 and 2. Table (3) shows the
$W_{m\mu}^{M}$ values considering this approach for the example in equation
(17).

| Milestone Years $\mu$
---|---
Milestone Year $m$ |  0 |  2 |  5
0 |  1 |  1/2 |  0
2 |  0 |  1/2 |  1/2
5 |  0 |  0 |  1/2
Table 3: Weight of investment milestone year $m$ on operational milestone year
$\mu$ ($W_{m\mu}^{M}$)

Another approach to determining the $W_{m\mu}^{M}$ values could be to consider
an efficiency on the capacity factor (or degradation) that changes over the
milestone years, such as the investments on milestone years with more
efficiency get a higher proportion of the operational production. So, for
instance, if years 2 and 5 are twice efficient as the technology invested in
year 0, then $W_{0,2}^{M}=1/3$ and $W_{2,2}^{M}=2/3$. The assumption here is
that all investments will be made. In other words, the distribution among the
investment decisions is predefined before the optimization decision. However,
if $x_{0}=0$ and $x_{2}\neq 0$, then we are underestimating the cost recovery
for year 2. The proposed formulation in Section 5.3 shows a workaround to
avoid this assumption.

For the parameter $W_{my\mu}^{Y}$, a simple approach distributes it linearly
between milestone years $m$ and $\mu$. Table (4) shows the values for the
example we have been considering before. For instance, the non-modelled year 1
is between milestone years 0 and 2; therefore, $W_{0,1,0}^{Y}=1/2$ and
$W_{0,1,2}^{Y}=1/2$, meaning that we represent the operational cost of year 1
with half the value in year 0 and half in year 2. Notice that the values for
$W_{0,2,2}^{Y}=W_{0,3,2}^{Y}=W_{0,4,2}^{Y}=1$ because the lifetime of the
investment on year 0 ends on year 4; so the operational production of
milestone year 5 does not consider the investment on year 0; see equation
(17).

|  | Milestone Year $\mu$
---|---|---
Milestone Year $m$ | Year $y$ |  0 |  2 |  5
0 | 0 |  1 |  0 |  0
| 1 |  1/2 |  1/2 |  0
| 2 |  0 |  1 |  0
| 3 |  0 |  1 |  0
| 4 |  0 |  1 |  0
| 5 |  0 |  0 |  0
2 | 2 |  0 |  1 |  0
| 3 |  0 |  2/3 |  1/3
| 4 |  0 |  1/3 |  2/3
| 5 |  0 |  0 |  1
5 | 5 |  0 |  0 |  1
Table 4: Weight of milestone year $m$ on year $y$ and milestone year $\mu$
($W_{my\mu}^{Y}$)

Using the information of Table 3, Table 4, $LT=5$, and the equation (16), the
total operational cost will be:

$\displaystyle C^{O}=$
$\displaystyle\sum_{m\in\mathcal{M}}\sum_{y=m}^{\min(m+LT-1,Y-1)}\frac{1}{(1+R)^{y}}\sum_{\mu\in\mathcal{M}}{\color[rgb]{0,0,1}W_{my\mu}^{Y}}\cdot
C_{\mu}^{op}\sum_{k}W_{\mu
k}^{op}\cdot\sum_{t}{\color[rgb]{.75,0,.25}W_{m\mu}^{M}}\cdot p_{\mu kt}$
$\displaystyle=C_{0}^{op}\sum_{k}W_{0,k}^{op}\cdot\sum_{t}p_{0,kt}$
$\displaystyle+\frac{1}{(1+R)^{1}}\frac{1}{2}C_{0}^{op}\sum_{k}W_{0,k}^{op}\cdot\sum_{t}p_{0,kt}+\frac{1}{(1+R)^{1}}\frac{1}{2}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{2,kt}$
$\displaystyle+\frac{1}{(1+R)^{2}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{2,kt}$
$\displaystyle+\frac{1}{(1+R)^{3}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{2,kt}$
$\displaystyle+\frac{1}{(1+R)^{4}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{2,kt}$
$\displaystyle+\frac{1}{(1+R)^{2}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{2,kt}$
$\displaystyle+\frac{1}{(1+R)^{3}}\frac{2}{3}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{2,kt}+\frac{1}{(1+R)^{3}}\frac{1}{3}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{5,kt}$
$\displaystyle+\frac{1}{(1+R)^{4}}\frac{1}{3}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{2,kt}+\frac{1}{(1+R)^{4}}\frac{2}{3}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{5,kt}$
$\displaystyle+\frac{1}{(1+R)^{5}}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{5,kt}$
$\displaystyle+\frac{1}{(1+R)^{5}}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}\frac{1}{2}p_{5,kt}$

Here the first five terms correspond to the operational cost of investment
made in year 0, the following four terms to the investment in milestone year
2, and the last one to milestone year 5. Therefore, it improves the previous
proposal in Section 5.1 by correctly accounting for each investment’s
operational cost in the milestone year, avoiding overestimating the
operational cost. However, as mentioned before, the central assumption in this
proposal is the predefined distribution of the operational costs before the
optimization, which could potentially lead to an underestimation of the
operational cost.

###  5.3 Proposal 3 - New operational variable

Although proposals in Sections 5.1 and 5.2 improve the standard formulation in
Section 4.2, they either overestimate or underestimate the operational costs
in the objective function. The main issue is that consolidating the production
data into a single variable per year poses a challenge in determining the
production outcome for each investment decision in advance. This is because
investment decisions are endogenous, resulting from the optimization problem.
Therefore, a more detailed formulation is required to accurately determine the
production of each investment decision; see equation (18).

$\displaystyle\min_{x_{m},p_{m\mu kt},p_{\mu kt}}\quad C^{I}+C^{O}$ (18)
$\displaystyle s.t.$ $\displaystyle p_{m\mu kt}\leq x_{m}\quad\forall
m\forall\mu\forall k\forall t|m\leq\mu\leq min(m+LT-1,Y-1)$ $\displaystyle
p_{\mu kt}=\sum_{m\geq\max(\mu-LT+1,0)}^{\mu}p_{m\mu kt}\quad\forall\mu\forall
k\forall t$

The new variable $p_{m\mu kt}$ represents the production imputed to the
investment in year $m$ for period $t$ and representative $k$ in operational
year $\mu$. In the example, it will be:

$\displaystyle p_{0,0,k,t}$ $\displaystyle\leq x_{0}\quad\forall k\forall t$
$\displaystyle p_{0,2,k,t}$ $\displaystyle\leq x_{0}\quad\forall k\forall t$
$\displaystyle p_{2,2,k,t}$ $\displaystyle\leq x_{2}\quad\forall k\forall t$
$\displaystyle p_{2,5,k,t}$ $\displaystyle\leq x_{2}\quad\forall k\forall t$
$\displaystyle p_{5,5,k,t}$ $\displaystyle\leq x_{5}\quad\forall k\forall t$
$\displaystyle p_{0,k,t}$ $\displaystyle=p_{0,0,k,t}\quad\forall k\forall t$
$\displaystyle p_{2,k,t}$ $\displaystyle=p_{0,2,k,t}+p_{2,2,k,t}\quad\forall
k\forall t$ $\displaystyle p_{5,k,t}$
$\displaystyle=p_{2,5,k,t}+p_{5,5,k,t}\quad\forall k\forall t$

Then the operational cost can be easily determined using the new variables and
the parameter $W_{my\mu}^{Y}$ as in equation (19).

$C^{O}=\sum_{m\in\mathcal{M}}\sum_{y=0}^{Y-1}\frac{1}{(1+R)^{y}}\sum_{\mu=m}^{min(m+LT-1,Y-1)}{\color[rgb]{0,0,1}W_{my\mu}^{Y}}\cdot
C_{m}^{op}\sum_{k}W_{mk}^{op}\cdot\sum_{t}p_{m\mu kt}$ (19)

For the example we have been using before, equation (19) becomes:

$\displaystyle C^{O}$
$\displaystyle=C_{0}^{op}\sum_{k}W_{0,k}^{op}\cdot\sum_{t}p_{0,0,k,t}$
$\displaystyle+\frac{1}{(1+R)^{1}}\frac{1}{2}C_{0}^{op}\sum_{k}W_{0,k}^{op}\cdot\sum_{t}p_{0,0,k,t}+\frac{1}{(1+R)^{1}}\frac{1}{2}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{0,2,k,t}$
$\displaystyle+\frac{1}{(1+R)^{2}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{0,2,k,t}$
$\displaystyle+\frac{1}{(1+R)^{3}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{0,2,k,t}$
$\displaystyle+\frac{1}{(1+R)^{4}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{0,2,k,t}$
$\displaystyle+\frac{1}{(1+R)^{2}}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{2,2,k,t}$
$\displaystyle+\frac{1}{(1+R)^{3}}\frac{2}{3}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{2,2,k,t}+\frac{1}{(1+R)^{3}}\frac{1}{3}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}p_{2,5,k,t}$
$\displaystyle+\frac{1}{(1+R)^{4}}\frac{1}{3}C_{2}^{op}\sum_{k}W_{2,k}^{op}\cdot\sum_{t}p_{2,2,k,t}+\frac{1}{(1+R)^{4}}\frac{2}{3}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}p_{2,5,k,t}$
$\displaystyle+\frac{1}{(1+R)^{5}}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}p_{2,5,k,t}$
$\displaystyle+\frac{1}{(1+R)^{5}}C_{5}^{op}\sum_{k}W_{5,k}^{op}\cdot\sum_{t}p_{5,5,k,t}$

Here it is possible to see that all the issues stated in Section 5 are solved.
In addition, this proposal avoids both overestimating and underestimating the
operational costs. However, it does all this at the cost of having more
variables in the optimization problem. For instance, there are five more
variables in the numerical example.

##  6 Formulations Summary

Throughout this document, we analyzed various formulations and modifications
to accurately account for the investments and operational costs of multi-year
investment optimization problems in energy systems. We have demonstrated that
both the total investment cost method and the annualized investment cost
method lead to the same representation if executed carefully. Additionally,
the operational cost representation presents challenges in avoiding issues
with the value of money over time and the asset’s lifetime. The proposals in
Section 5 provide a comprehensive representation that improves the current
standard methods in the literature. Here is a summary of all the discussed
formulations:

  * •

Weights on milestone years (standard approach)

    * –

Equations: (12) and (13)

    * –

PROS: Simple and straightforward to use.

    * –

CONS: It does not account appropriately for the investment costs in the
annualized formulation nor the value of money over time

  * •

Proposal 1 - Mapping investment milestone years to all years ($W_{my}^{O}$)

    * –

Equations: (14) and (15)

    * –

PROS: It accounts properly for the investment costs in the annualized
formulation as well as the value of money over time

    * –

CONS: It overestimates the operational costs

  * •

Proposal 2 - Mapping investment milestone years to all years and operational
milestone years ($W_{my\mu}^{Y}$ and $W_{m\mu}^{M}$)

    * –

Equations: (14) and (15), but changing the $C^{O}$ for the one in (16)

    * –

PROS: It accounts properly for the investment costs in the annualized
formulation as well as the value of money over time

    * –

CONS: It underestimates the operational costs

  * •

Proposal 3 - New operational variable

    * –

Equations: (14) and (15), but changing the $C^{O}$ for the one in (19)

    * –

PROS: It accounts for all the costs in the objective function properly

    * –

CONS: More variables in the optimization problem

###  6.1 Final Comments

In the event that the investment is a fixed value and the optimization problem
is centred solely on operational costs, Proposal 2 can be executed without
underestimations on production costs.

By extending the values of weights beyond the last year of modelling in any of
the proposals, we can consider new production costs instead of solely relying
on the recovery assumption. This approach broadens the scope of economic
assumptions in analysing assets barely produced in the last modelled year.
Furthermore, it enables us to have a lengthier cash flow projection, which is
instrumental in double-checking the profitability of the investment without
necessarily increasing the number of milestone years.

Finally, this document’s future versions will analyse the formulation changes
when considering decommissioning variables and vintage investment
representation.

## References

  * [1] Efraim Centeno, Sonja Wogrin and Diego Tejada-Arango  “On the proper use of net present value for capacity expansion planning with renewable generation”, Working paper IIT-18-128A, 2018
  * [2] T Brown  “Multi-Horizon Planning with Perfect Foresight”, Preprint submitted to Elsevier, 2020  URL: https://nworbmot.org/energy/multihorizon.pdf