# Tunable and efficient ultraviolet generation in nanophotonic lithium niobate

Emily Hwang  1,$\dagger$ Nathan Harper  2,$\dagger$ Ryoto Sekine  3 Luis
Ledezma  3,4 Alireza Marandi  1,3 and Scott K. Cushing2  1Department of
Applied Physics and Materials Science, California Institute of Technology,
Pasadena, California, 91125, USA
2Division of Chemistry and Chemical Engineering, California Institute of
Technology, Pasadena, California 91125, USA
3Department of Electrical Engineering, California Institute of Technology,
Pasadena, California 91125, USA
4 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA
91109, USA
$\dagger$ These authors contributed equally to this work.  *
<EMAIL_ADDRESS>
###### Abstract

On-chip ultraviolet sources are of great interest for building compact and
scalable atomic clocks, quantum computers, and spectrometers; however, few
material platforms are suitable for integrated ultraviolet light generation.
Of these materials, thin-film lithium niobate is the most competitive due to
its ability to be quasi-phase-matched, optical confinement, and nonlinear
properties. Here, we present efficient (197 $\pm$ 5 %/W/cm2) second harmonic
generation of UV-A light in a periodically poled lithium niobate nanophotonic
waveguide. We achieve on-chip ultraviolet powers of 30 $\mu$W, demonstrating
the potential for compact frequency-doubling of common near-IR laser diodes.
By using a large cross section waveguide (600 nm film thickness), we achieve
insensitivity to fabrication errors, and can attain first-order quasi-phase-
matching with relatively long poling periods (>1.5 $\mu$m). The device also
demonstrates linear wavelength tunability using temperature. By varying the
poling period, we have achieved the shortest reported wavelength (355 nm)
generated through frequency doubling in thin-film lithium niobate. Our results
open up new avenues to realize ultraviolet on-chip sources and chip-scale
photonics.

††journal: opticajournal††articletype: Letter

The field of integrated nonlinear optics has grown dramatically during the
past decade due to the development and commercial availability of thin-film
lithium niobate[1, 2]. In passive nonlinear devices, thin-film lithium niobate
(TFLN) excels in efficient frequency conversion and quantum state generation
from the visible to the infrared (IR)[3, 4]. The strong mode confinement of
single-pass, low-loss[5] nanophotonic waveguides and quasi-phase matched
interactions utilizing lithium niobate’s largest second-order nonlinear
optical tensor element have resulted in record-breaking efficiencies in
applications such as second harmonic generation (SHG)[6, 7], supercontinuum
generation[8], difference frequency generation[9], parametric
amplification[10], and parametric downconversion[11]. Similarly, active
devices such as modulators[12], electro-optic frequency combs[13], and
femtosecond pulse generators[14] show impressive performance in compact form-
factors due to lithium niobate’s large electro-optic tensor elements. However,
there is still significant room for lithium niobate’s use in ultraviolet (UV)
photonics[15], with applications such as UV-visible spectroscopy,
optogenetics, high-resolution microscopy, UV security banknote features, laser
cooling[16], atomic clocks[17, 18], and quantum computing[19].

Although lithium niobate has been extensively studied in the IR, and
comparatively less so in the visible, UV devices have remained rare to date.
The few reported lithium niobate devices for UV generation have been limited
to metasurfaces[20], nanoparticles[21], and large micromachined or channel
waveguides[22, 23], and therefore do not take advantage of the high mode
confinement and efficiency of lithium niobate in a nanophotonic platform. TFLN
has yet to be well studied in the UV due to the ultra-short poling periods
required to overcome the high dispersion in waveguides at short wavelengths,
as well as material and scattering loss. Although lithium niobate’s band gap
is near 315 nm (3.95 eV), an exponentially decaying Urbach absorption tail
persists towards the visible due to defects in the crystal structure[24], and
impurity ion (Cu+, Fe2+) resonances can cause additional loss[25, 26].
Additionally, losses at the waveguide sidewalls increase at shorter
wavelengths due to surface imperfection Rayleigh scattering, which scales as
$\lambda^{-4}$[27]. In spite of these difficulties, there is much to gain by
extending the spectral coverage of TFLN frequency conversion to the UV.
Notably, near-IR laser diodes, which can be frequency doubled, are
considerably more accessible than UV laser diodes and gas lasers[28]. Among
other nanophotonic material platforms, only aluminum nitride (AlN) has been
significantly investigated for waveguided second harmonic UV generation;
however, AlN lacks ferroelectricity and is therefore incapable of periodic
poling. The lateral polar structures used to achieve quasi-phase matching in
AlN are highly scattering, resulting in much lower conversion efficiencies
(<1%)[29] compared to lithium niobate devices. Other potential UV platforms
(lithium tantalate[30], BBO[31], LBGO[32]) have yet to be thoroughly explored
in a thin film nanophotonic platform. Lithium niobate remains superior to
these materials with its combination of low-loss waveguides, high second-order
nonlinear response, ferroelectric poling for quasi-phase matching, and, for
the case of thin-film lithium tantalate, commercial accessibility[33].

Here, we produce 30 $\mu$W of efficient (197 $\pm$ 5 %/W/cm2) second harmonic
generation of UV light (386.5 nm) with periodically poled lithium niobate
(PPLN) rib waveguides (Fig. 1a). In addition to being among the first reported
UV SHG devices in thin-film lithium niobate, the waveguides are designed to be
robust against fabrication errors. The devices exhibits wavelength tunability
through temperature and poling period, and is capable of UV SHG at the lowest
wavelengths tested (710 nm frequency-doubled to 355 nm).

The optimal waveguide geometry was determined by minimizing the phase matching
sensitivity to thickness variations in the LN thin film while maintaining the
SHG bandwidth at the center wavelength of the laser used in this paper (773
nm). To date, variations in the thin film thickness of even 1 Å are sufficient
to disrupt phase matching in visible SHG, limiting the effective interaction
length and chip-to-chip repeatability[1, 34]. To quantify the phase matching
sensitivity, the field profiles and effective indices of the guided modes were
simulated (Lumerical MODE) using the bulk Sellmeier coefficients of lithium
niobate[35] and SiO2[36] with the geometric parameters shown in Fig. 1b and a
sidewall angle of 60°, which is consistent with the fabrication process. Only
the fundamental quasi-transverse electric (TE) modes of the first harmonic
(FH) and second harmonic (SH) (Fig. 1b) were considered since they access
lithium niobate’s largest nonlinear tensor element (d33 = 25 pm/V)[37]. The
first-order wavelength sensitivity and SHG bandwidth were calculated using a
numerical first derivative with respect to film thickness over a range of
waveguide geometries for 200 nm and 600 nm thicknesses (Fig. 1c). For a
waveguide top width of 1.5 $\mu$m and etch depth of 375 nm (Fig. 1d-e) with a
600 nm film, the wavelength sensitivity to film thickness is
$\frac{d\lambda}{dT}=$ 0.08 nm/nm with a 3.4 pm SHG bandwidth, and the
effective refractive indices ($n_{\text{eff,SH}}$ = 2.32, $n_{\text{eff,FH}}$
= 2.10) result in a quasi-phase matching poling period of
$\Lambda=\frac{\lambda_{\text{SH}}}{(\Delta n_{\text{eff}})}=1.8\mu\text{m}$
(Fig. 1f). The larger cross section of the 600 nm film thickness allows for
this relatively long poling period, in comparison to a 1.1 $\mu$m period for a
200 nm film with the same aspect ratio and wavelength. By using a waveguide
with a larger cross section, we achieve insensitivity to the film thickness
and fabrication errors without sacrificing the phase matching bandwidth.

Figure 1: (a) Schematic of the PPLN waveguide. (b) Mode profiles of the
fundamental TE mode at the first and second harmonics. (c) Sensitivity of the
SHG center wavelength to the thin film thickness as a function of the
waveguide geometry for a 600 nm (top) and a 200 nm thin film thickness
(bottom); local areas that minimize the sensitivity also decrease the phase
matching bandwidth. Two dimensional sweeps of the (d) film thickness and etch
depth and (e) top width and etch depth to vary the SHG center wavelength
(contour lines). (f) Poling period and effective refractive indices of the
first and second harmonic fundamental TE modes as a function of wavelength.

The devices were fabricated from a 5% MgO-doped X-cut thin-film lithium
niobate on insulator wafer (NANOLN), which consists of 600 nm of lithium
niobate bonded to 2 um of silicon dioxide on a 0.4 mm silicon substrate.
Periodically poled waveguides with a 1.5 $\mu$m top width and 375 nm etch
depth were fabricated following Ref. [10]. Each waveguide had a 7 mm poled
length, with poling periods ranging from 1.55 to 2 $\mu$m. The waveguide etch
depth and sidewall angle were verified through atomic force microscopy, and
the formation of poled domains was measured with second harmonic
microscopy[38].

The SHG from the PPLN devices was characterized by the optical setup in Fig.
2a. The output from a tunable continuous-wave (CW) single-frequency laser
(Velocity TLD-6712, 765-781 nm) passed through an optical isolator and a
variable neutral density filter to adjust the input power. An achromatic half-
wave plate (Thorlabs AHWP10M-980) aligned the input polarization to the
optical axis of the chip to maximize SHG power. The first harmonic was coupled
to the waveguide using an AR-coated aspheric lens (focal length 1.5 mm,
Thorlabs C140TMD-B). The waveguide output was collimated by another aspheric
lens (focal length 1.5 mm, Thorlabs C140TMD-A) and collected by a high-OH
multimode fiber (Thorlabs M122L01). The second harmonic was monitored using an
optical spectrum analyzer (Yokogawa AQ6374) with a passband bandwidth of 5 nm
around the second harmonic center wavelength to remove residual first
harmonic.

Figure 2:  (a) Optical setup for on-chip SHG characterization. ND, variable
neutral density filter; HWP, half wave plate; L, aspheric lens; MMF, multi-
mode fiber; OSA, optical spectrum analyzer. (b) On-chip SHG output power
measured as a function of the on-chip input power and experimental fit to
calculate the efficiency. (c) Temperature tuning of the SHG center wavelength,
measured as the weighted average of the spectrum with corresponding standard
deviations and experimental fit. (d) Normalized spectra of the poling period
sweep.

The off-chip power of the first harmonic was varied from 6 to 30 mW at 773.1
nm, and the resulting on-chip SHG power is plotted in Fig. 2b. We achieved a
maximum on-chip UV power of 30 $\mu$W, and the slope of the least squares fit
of the SHG output gives an efficiency of 197 $\pm$ 5 %/W/cm2. The power law
scaling suggests that photorefraction and pump depletion are not significant
at the maximum power tested. The calculated efficiency and on-chip powers take
the coupling losses into account, which were 22% transmission per facet for
the first harmonic and 10% transmission per facet for the second harmonic.
These values were determined by measuring the overall transmission at 780 nm
and 405 nm, and assuming that the input and output coupling losses are equal.

The temperature and poling period dependence of the phase matched wavelength
is plotted in Fig. 2c and 2d, respectively. To measure the temperature
dependence, a nonlinear crystal oven (HC Photonics TC038-PC) heated the PPLN
waveguides from 25 to 100 °C. The least squares slope (34 $\pm$ 2 pm/°C) is in
close agreement with the theoretically calculated value of 32 pm/°C. By
varying the poling period, we also obtained SHG spectra extending down to 355
nm. A tunable CW Ti:Sapphire oscillator (Spectra-Physics Tsunami, 700-1100 nm)
was used in a similar scheme to Fig. 2a to measure the SHG from several PPLN
waveguides with poling periods as short as 1.35 $\mu$m (Fig. 2d). The fitted
slope of the SHG center wavelength with the poling period (79.6 $\pm$ 0.2
nm/$\mu$m) matches the theoretical slope of 81.1 nm/$\mu$m. This agreement in
both the temperature and poling period demonstrate that the given Sellmeier
coefficients[35] predict the temperature dependence and group velocity
mismatch of lithium niobate in the UV relatively well, despite the limited
refractive index data at these wavelengths. Unlike the initial measurements
using the single-frequency laser, accurate efficiency data could not be
extracted because the Ti:Sapphire oscillator linewidth (0.3 nm) is orders of
magnitude larger than the phase matching bandwidth (3.4 pm). However, these
spectra demonstrate that SHG is possible even closer to the 315 nm band gap of
lithium niobate than what has previously been demonstrated, and that the thin
film lithium niobate platform is able to phase match the full range of a
Ti:Sapphire laser.

Although the SHG temperature and poling period tuning agree well with theory,
the resulting 197 $\pm$ 5 %/W/cm2 efficiency measured from the single-
frequency laser at 773.1 nm is lower than the calculated theoretical
efficiency[7] of 18,100 %/W/cm2. The experimental SHG spectrum (Fig. 3a) was
also measured by sweeping the first harmonic wavelength from 771 to 775 nm
with a constant input power of 30 mW. The spectrum exhibits multiple peaks
over a 0.5 nm bandwidth, deviating from the theoretical sinc2 line shape and
3.4 pm FWHM. Although the waveguides in this work support multiple modes at
both the first and second harmonic, significant phase matching between higher-
order modes is unlikely due to large momentum mismatches or poor modal
overlap. Furthermore, the shape of the transfer function does not change when
the second harmonic is collected with a single-mode fiber (Thorlabs
P1-305A-FC-1), which is expected to occur if higher-order second harmonic
modes were present.

The discrepancies in the efficiency, bandwidth, and spectral shape are likely
caused by index variations, potentially from thickness variations, thermal
gradients, induced absorption, or photorefraction. Index variations preserve
the area under the curve of the SHG transfer function, which allows us to
determine the contribution of the index variations to the experimentally lower
efficiency. An integral of the spectral efficiency yields an area of 31.8
%/W/cm${}^{2}\cdot$nm, and the transfer function with the same area and no
index variations has a peak efficiency of 861 %/W/cm2. The poling quality of
this device, which exhibited significant domain widening due to the short
poling period, also contributes to the lowered efficiency. The duty cycle is
estimated to be 90% from second harmonic microscopy images, which lowers the
theoretical efficiency to 1810 %/W/cm2. The remaining discrepancies can be
explained by loss at the second harmonic and asymmetry in the input and output
coupling efficiency.

Figure 3:  (a) Experimental and index variation-fitted SHG spectra. (b) Phase
mismatch error represented by Eq. 2 and corresponding thickness error against
propagation distance.

The magnitude of the index variations through the waveguide are estimated
through simulation. Without any knowledge of the SHG phase relative to the
first harmonic, the magnitude or position of the index variations cannot be
directly calculated. However, the SHG spectrum can be approximated by the
following expression for $\eta(\omega)$[39]:

$\eta(\omega)=\bigg{|}\int_{0}^{L}\text{exp}\left(i\int_{0}^{z}(\delta\beta(\xi)+\delta\beta(\omega))\,d\xi\right)\,dz\bigg{|}^{2}$
(1)

Where the phase mismatch error $\delta\beta(z)$ is expressed along the length
of the waveguide as a cubic polynomial. The coefficients of $\delta\beta(z)$
are estimated by minimizing the squared residuals between the predicted and
experimental spectra using a particle swarm optimization algorithm followed by
a gradient descent[40, 41, 42], which yields the following expression for
$\delta\beta(z)$ (Fig. 3b):

$\delta\beta(z)=0.0636z^{3}-0.578z^{2}+0.544z+3.66$ (2)

The blue trace in Fig. 3a is the predicted spectrum using Eq. 2. A thickness
error of $\pm$5 nm throughout the waveguide explains the full spread in the
experimental spectrum if the index variations are attributed only to
thickness. A film thickness discrepancy of this magnitude is reasonable given
the height measurements performed by NANOLN; however, this error is an upper
estimate given that the index variations could be a combination of the
thickness, photorefraction, UV-induced infrared absorption[43], and
temperature gradients. Although this thickness variation could be the cause of
the lowered experimental efficiency, the same amount of thickness variation in
a smaller film thickness could potentially disrupt phase matching entirely.
The large cross sections of the waveguides makes the devices robust to the
total film thickness, and therefore allows us to still achieve phase matching
despite the lowered efficiency. Through this simple approximation, we believe
that future devices can be optimized by estimating deviations from the ideal
index, and that this method can serve as a means of monitoring
photorefraction.

In conclusion, we have demonstrated fabrication-insensitive temperature-
tunable ultraviolet light generation in an integrated thin-film lithium
niobate waveguide. We have observed an experimental SHG efficiency of 197
$\pm$ 5 %/W/cm2, with discrepancies from the theoretical efficiency that can
be explained by the poling duty cycle and index variations in the waveguides.
As of this publication, this is the shortest wavelength (355-386 nm) produced
through second harmonic generation with periodically poled thin-film lithium
niobate waveguides. Our work opens up opportunities to realize efficient
frequency-doubled chip-scale ultraviolet laser diodes for UV integrated
photonics, with applications spanning spectroscopy, atomic physics, and
quantum science.

Funding Department of Energy (DE-SC0020151); KNI-Wheatley Scholar in
Nanoscience; National Science Foundation (DGE-1745301); Rothenberg Innovation
Initiative.

Acknowledgments The device fabrication was performed at the Kavli Nanoscience
Institute at Caltech. The authors thank Robert Gray for his experimental
support. E.H. was supported by the National Science Foundation Graduate
Research Fellowship Program under Grant no. DGE-1745301. Any opinion,
findings, and conclusions or recommendations expressed in this material are
those of the authors and do not necessarily reflect the views of the National
Science Foundation. N.H. was supported by the Department of Defense (DoD)
through the National Defense Science and Engineering Graduate (NDSEG)
Fellowship Program.

Disclosures The authors declare no conflicts of interest.

Data Availability Statement Data underlying the results presented in this
paper are not publicly available at this time but may be obtained from the
authors upon reasonable request.

## References

  * [1] A. Boes, L. Chang, C. Langrock, M. Yu, M. Zhang, Q. Lin, M. Lončar, M. Fejer, J. Bowers, and A. Mitchell, “Lithium niobate photonics: Unlocking the electromagnetic spectrum,” Science 379, eabj4396 (2023).
  * [2] D. Zhu, L. Shao, M. Yu, R. Cheng, B. Desiatov, C. J. Xin, Y. Hu, J. Holzgrafe, S. Ghosh, A. Shams-Ansari, E. Puma, N. Sinclair, C. Reimer, M. Zhang, and M. Lončar, “Integrated photonics on thin-film lithium niobate,” Adv. Opt. Photon. 13, 242 (2021).
  * [3] R. Nehra, R. Sekine, L. Ledezma, Q. Guo, R. M. Gray, A. Roy, and A. Marandi, “Few-cycle vacuum squeezing in nanophotonics,” Science 377, 1333–1337 (2022).
  * [4] M. Leidinger, S. Fieberg, N. Waasem, F. Kühnemann, K. Buse, and I. Breunig, “Comparative study on three highly sensitive absorption measurement techniques characterizing lithium niobate over its entire transparent spectral range,” Opt. Express 23, 21690 (2015).
  * [5] A. Shams-Ansari, G. Huang, L. He, Z. Li, J. Holzgrafe, M. Jankowski, M. Churaev, P. Kharel, R. Cheng, D. Zhu, N. Sinclair, B. Desiatov, M. Zhang, T. J. Kippenberg, and M. Lončar, “Reduced material loss in thin-film lithium niobate waveguides,” APL Photonics 7, 081301 (2022).
  * [6] T. Park, H. S. Stokowski, V. Ansari, T. P. McKenna, A. Y. Hwang, M. M. Fejer, and A. H. Safavi-Naeini, “High-efficiency second harmonic generation of blue light on thin-film lithium niobate,” Opt. Lett. 47, 2706 (2022).
  * [7] M. Jankowski, C. Langrock, B. Desiatov, A. Marandi, C. Wang, M. Zhang, C. R. Phillips, M. Lončar, and M. M. Fejer, “Ultrabroadband nonlinear optics in nanophotonic periodically poled lithium niobate waveguides,” Optica 7, 40 (2020).
  * [8] C. Wang, C. Langrock, A. Marandi, M. Jankowski, M. Zhang, B. Desiatov, M. M. Fejer, and M. Lončar, “Ultrahigh-efficiency wavelength conversion in nanophotonic periodically poled lithium niobate waveguides,” Optica 5, 1438 (2018).
  * [9] J. Mishra, T. P. McKenna, E. Ng, H. S. Stokowski, M. Jankowski, C. Langrock, D. Heydari, H. Mabuchi, M. M. Fejer, and A. H. Safavi-Naeini, “Mid-infrared nonlinear optics in thin-film lithium niobate on sapphire,” Optica 8, 921 (2021).
  * [10] L. Ledezma, R. Sekine, Q. Guo, R. Nehra, S. Jahani, and A. Marandi, “Intense optical parametric amplification in dispersion-engineered nanophotonic lithium niobate waveguides,” Optica 9, 303 (2022).
  * [11] J. Zhao, C. Ma, M. Rüsing, and S. Mookherjea, “High Quality Entangled Photon Pair Generation in Periodically Poled Thin-Film Lithium Niobate Waveguides,” Phys. Rev. Lett. 124, 163603 (2020).
  * [12] C. Wang, M. Zhang, X. Chen, M. Bertrand, A. Shams-Ansari, S. Chandrasekhar, P. Winzer, and M. Lončar, “Integrated lithium niobate electro-optic modulators operating at CMOS-compatible voltages,” Nature 562, 101–104 (2018).
  * [13] M. Zhang, B. Buscaino, C. Wang, A. Shams-Ansari, C. Reimer, R. Zhu, J. M. Kahn, and M. Lončar, “Broadband electro-optic frequency comb generation in a lithium niobate microring resonator,” Nature 568, 373–377 (2019).
  * [14] M. Yu, D. Barton III, R. Cheng, C. Reimer, P. Kharel, L. He, L. Shao, D. Zhu, Y. Hu, H. R. Grant, L. Johansson, Y. Okawachi, A. L. Gaeta, M. Zhang, and M. Lončar, “Integrated femtosecond pulse generator on thin-film lithium niobate,” Nature 612, 252–258 (2022).
  * [15] M. Soltani, R. Soref, T. Palacios, and D. Englund, “AlGaN/AlN integrated photonics platform for the ultraviolet and visible spectral range,” Opt. Express 24, 25415 (2016).
  * [16] K. Toyoda, A. Miura, S. Urabe, K. Hayasaka, and M. Watanabe, “Laser cooling of calcium ions by use of ultraviolet laser diodes: significant induction of electron-shelving transitions,” Opt. Lett. 26, 1897 (2001).
  * [17] K. Sugiyama and J. Yoda, “Production of YbH${}^{\textrm{+}}$ by chemical reaction of Yb${}^{\textrm{+}}$ in excited states with H${}_{\textrm{2}}$ gas,” Phys. Rev. A 55, R10–R13 (1997).
  * [18] A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. Schmidt, “Optical atomic clocks,” Rev. Mod. Phys. 87, 637–701 (2015).
  * [19] M. Saffman, “Quantum computing with atomic qubits and Rydberg interactions: progress and challenges,” J. Phys. B: At. Mol. Opt. Phys. 49, 202001 (2016).
  * [20] J. Ma, F. Xie, W. Chen, J. Chen, W. Wu, W. Liu, Y. Chen, W. Cai, M. Ren, and J. Xu, “Nonlinear Lithium Niobate Metasurfaces for Second Harmonic Generation,” Laser & Photonics Reviews 15, 2000521 (2021).
  * [21] F. Timpu, J. Sendra, C. Renaut, L. Lang, M. Timofeeva, M. T. Buscaglia, V. Buscaglia, and R. Grange, “Lithium Niobate Nanocubes as Linear and Nonlinear Ultraviolet Mie Resonators,” ACS Photonics 6, 545–552 (2019).
  * [22] J. Rutledge, A. Catanese, D. D. Hickstein, S. A. Diddams, T. K. Allison, and A. S. Kowligy, “Broadband ultraviolet-visible frequency combs from cascaded high-harmonic generation in quasi-phase-matched waveguides,” J. Opt. Soc. Am. B 38, 2252 (2021).
  * [23] T. Sugita, K. Mizuuchi, Y. Kitaoka, and K. Yamamoto, “Ultraviolet Light Generation in a Periodically Poled MgO:LiNbO${}_{\textrm{3}}$ Waveguide,” Jpn. J. Appl. Phys. 40, 1751 (2001).
  * [24] R. Bhatt, I. Bhaumik, S. Ganesamoorthy, A. K. Karnal, M. K. Swami, H. S. Patel, and P. K. Gupta, “Urbach tail and bandgap analysis in near stoichiometric LiNbO${}_{\textrm{3}}$ crystals,” Phys. Status Solidi A 209, 176–180 (2012).
  * [25] J. R. Schwesyg, M. C. C. Kajiyama, M. Falk, D. H. Jundt, K. Buse, and M. M. Fejer, “Light absorption in undoped congruent and magnesium-doped lithium niobate crystals in the visible wavelength range,” Appl. Phys. B 100, 109–115 (2010).
  * [26] M. V. Ciampolillo, A. Zaltron, M. Bazzan, N. Argiolas, and C. Sada, “Quantification of Iron (Fe) in Lithium Niobate by Optical Absorption,” Appl Spectrosc 65, 216–220 (2011).
  * [27] R. Maurer, “Glass fibers for optical communications,” Proc. IEEE 61, 452–462 (1973).
  * [28] S. M. N. Hasan, W. You, M. S. I. Sumon, and S. Arafin, “Recent Progress of Electrically Pumped AlGaN Diode Lasers in the UV-B and -C Bands,” Photonics 8, 267 (2021).
  * [29] D. Alden, T. Troha, R. Kirste, S. Mita, Q. Guo, A. Hoffmann, M. Zgonik, R. Collazo, and Z. Sitar, “Quasi-phase-matched second harmonic generation of UV light using AlN waveguides,” Appl. Phys. Lett. 114, 103504 (2019).
  * [30] J.-P. Meyn, C. Laue, R. Knappe, R. Wallenstein, and M. Fejer, “Fabrication of periodically poled lithium tantalate for UV generation with diode lasers,” Appl Phys B 73, 111–114 (2001).
  * [31] K. Devi, S. Parsa, and M. Ebrahim-Zadeh, “Continuous-wave, single-pass, single-frequency second-harmonic-generation at 266 nm based on birefringent-multicrystal scheme,” Opt. Express 24, 8763 (2016).
  * [32] N. Umemura, J. Hirohashi, Y. Nakahara, H. Oda, and Y. Furukawa, “Temperature-dependent quasi phase-matching properties of periodically poled LaBGeO ${}_{\textrm{5}}$,” Opt. Mater. Express 9, 2159 (2019).
  * [33] D. J. Blumenthal, “Photonic integration for UV to IR applications,” APL Photonics 5, 020903 (2020).
  * [34] A. A. Sayem, Y. Wang, J. Lu, X. Liu, A. W. Bruch, and H. X. Tang, “Efficient and tunable blue light generation using lithium niobate nonlinear photonics,” Appl. Phys. Lett. 119, 231104 (2021).
  * [35] O. Gayer, Z. Sacks, E. Galun, and A. Arie, “Temperature and wavelength dependent refractive index equations for MgO-doped congruent and stoichiometric LiNbO${}_{\textrm{3}}$,” Appl. Phys. B 91, 343–348 (2008).
  * [36] T. Toyoda and M. Yabe, “The temperature dependence of the refractive indices of fused silica and crystal quartz,” J. Phys. D: Appl. Phys. 16, L97–L100 (1983).
  * [37] I. Shoji, T. Kondo, A. Kitamoto, M. Shirane, and R. Ito, “Absolute scale of second-order nonlinear-optical coefficients,” J. Opt. Soc. Am. B 14, 2268 (1997).
  * [38] M. Rüsing, J. Zhao, and S. Mookherjea, “Second harmonic microscopy of poled x-cut thin film lithium niobate: Understanding the contrast mechanism,” Journal of Applied Physics 126, 114105 (2019).
  * [39] S. Helmfrid and G. Arvidsson, “Influence of randomly varying domain lengths and nonuniform effective index on second-harmonic generation in quasi-phase-matching waveguides,” J. Opt. Soc. Am. B 8, 797 (1991).
  * [40] J. Kennedy and R. Eberhart, “Particle Swarm Optimization,” in _Proceedings of ICNN’95 - International Conference on Neural Networks,_ vol. 4 (Perth, Australia, 1995), pp. 1942–1948.
  * [41] E. Mezura-Montes and C. A. Coello Coello, “Constraint-handling in nature-inspired numerical optimization: Past, present and future,” Swarm and Evolutionary Computation 1, 173–194 (2011).
  * [42] M. E. H. Pedersen, “Good Parameters for Particle Swarm Optimization,” Tech. rep., Hvass Laboratories, Luxemborg (2010).
  * [43] A. M. Ali and Y. Huang, “Observation of ultra-long lifetime UV-induced infrared absorption in undoped and magnesium-doped congruent lithium niobate crystals,” Applied Surface Science 257, 5957–5960 (2011).


If the condition in Line 13 is satisfied in the $k$th iteration in the for-
loop of Lines 8-25, the necessary and sufficient conditions (i) and (ii) of
Case 1 will be satisfied, which implies that the LNP price $q^{l}$ lies in the
threshold price interval $(q_{a}^{k-1},q_{a}^{k})$. To further satisfy the
necessary and sufficient condition (iii) of Case 1, we execute Subroutine 3.1
in Line 14 and compute the LNP price $q^{l}$ that makes
$\tilde{R}_{a}^{pf}(q^{l})=0$ in Line 3 of Subroutine 3.1, where we utilize
the linear structure over the price interval as shown in Figure 9.

  * •

If the condition in Line 17 is satisfied in the $k$th iteration in the for-
loop of Lines 8-25, the necessary and sufficient conditions (i) and (ii) of
Case 2 will be both satisfied for the threshold price $q_{a}^{k}$. Thus, we
attain the LNP price $q^{l}$ at the threshold price $q_{a}^{k}$ as in Case 2.

  * •

If the condition in Line 21 is satisfied in the $k$th iteration in the for-
loop of Lines 8-25, the necessary and sufficient conditions (i) and (ii) of
Case 3 will be both satisfied for the threshold price $q_{a}^{k}$. Thus, we
attain the LNP price $q^{l}$ at the threshold price $q_{a}^{k}$ as in Case 3.

#### 2) The near-LNP price $\hat{q}^{l}$ under a positive $\varepsilon$

Since the aggregator cannot choose $\varepsilon=0$ due to the multi-optima
problem as we explained in Section IV.A, we cannot directly use the LNP price
$q^{l}$ in the limiting case where $\varepsilon$ approaches zero. Thus, we
propose an iterative procedure in Subroutines 3.1-3.3 of Algorithm 2 (for
Cases 1-3, respectively) to compute a near-LNP price $\hat{q}^{l}$ by
adjusting $\varepsilon$, such that i) we can still utilize the piecewise
linear solution structure in the limiting case; ii) we can compute a near-LNP
price to approximate the LNP price in the limiting case within a given
arbitrary accuracy. Furthermore, after computing the near-LNP price
$\hat{q}^{l}$ and $\varepsilon$, we can compute the investment decision
accordingly by Subroutine 1.

  * •

For Case 1, in Subroutine 3.1, we compute a small $\varepsilon$ in an
iterative fashion (by decreasing $\varepsilon$ iteratively), such that
${R_{a}^{pf}({q}^{l},\varepsilon)}$ is close to zero within an acceptable
accuracy $err_{4}$ (in Lines 4-8 of Subroutine 3.1). In this way, we use the
LNP price $q^{l}$ in the limiting case as the near-LNP price $\hat{q}^{l}$.

  * •

For Case 2, in Subroutine 3.2, the LNP price $q^{l}$ is attained at one
threshold price $q_{a}^{k}$ and $\tilde{R}_{a}^{pf}(q_{a}^{k-})=0$. Note that
we cannot directly choose the threshold price due to the discontinuity of
$\tilde{R}_{a}^{pf}(q)$. Instead, in Subroutine 3.2, since profit
$\tilde{R}_{a}^{pf}(q)$ is linearly increasing over $(q_{a}^{k-1},q_{a}^{k})$,
we let the aggregator choose a near-LNP price $\hat{q}^{l}$ slightly lower
than $q_{a}^{k}$, such that the profit ${R}_{a}^{pf}(\hat{q}^{l})$
approximates the profit $\tilde{R}_{a}^{pf}(q_{a}^{k-})$ within a given
accuracy $err_{3}$ (in Lines 3-4 of Subroutine 3.2).202020Recall that
${R}_{a}^{pf}(q)$ is equivalent to $\tilde{R}_{a}^{pf}(q)$ at non-threshold
prices. Then, the aggregator can compute a small $\varepsilon$ in an iterative
fashion so that the profit ${R}_{a}^{pf}(\hat{q}^{l},\varepsilon)$
approximates ${R}_{a}^{pf}(\hat{q}^{l})$ within the accuracy $err_{4}$ (in
Lines 5-9 of Subroutine 3.2).

  * •

For Case 3, in Subroutine 3.3, the LNP price $q^{l}$ is attained at one
threshold price $q_{a}^{k}$ and $\tilde{R}_{a}^{pf}(q_{a}^{k+})\geq 0$.
Similar to Case 2, we cannot directly choose the threshold price due to the
discontinuity of $\tilde{R}_{a}^{pf}(q)$. In Subroutine 3.3, since profit
$\tilde{R}_{a}^{pf}(q)$ is linearly increasing over $(q_{a}^{k},q_{a}^{k+1})$,
we let the aggregator choose a near-LNP price $\hat{q}^{l}$ slightly higher
than $q_{a}^{k}$, such that the profit ${R}_{a}^{pf}(\hat{q}^{l})$
approximates the profit $\tilde{R}_{a}^{pf}(q_{a}^{k+})$ within a given
accuracy $err_{3}$ (in Lines 3-4 of Subroutine 3.3). Then, the aggregator can
compute a small $\varepsilon$ in an iterative fashion so that the profit
${R}_{a}^{pf}(\hat{q}^{l},\varepsilon)$ approximates
${R}_{a}^{pf}(\hat{q}^{l})$ within the accuracy $err_{4}$ (in Lines 5-9 of
Subroutine 3.3).

Algorithm 2: Search of the near-LNP price $\hat{q}^{l}$

1: initialization: set iteration index $k=0$;

2: for each user $i\in\mathcal{I}$ in parallel do

3: Compute the set $\mathcal{Q}_{i}^{\omega}$ and $\mathcal{X}_{i}^{\omega}$
by Algorithm 3 of Appendix D;

4: Compute the limiting optimal charge/discharge decision corresponding to
each element in $\mathcal{X}_{i}^{\omega}$ by Algorithm 4 of Appendix H, for
all $\omega\in{\Omega}$;

5: Report all the computation results to the aggregator;

6: end for

7: The aggregator sorts the threshold price set
$\mathcal{Q}_{a}=\bigcup_{i,\omega}\mathcal{Q}_{i}^{\omega}$ by an increasing
order: $0=q_{a}^{0}<q_{a}^{1}<\cdots<q_{a}^{K_{a}}$.

8: for $k=0$ to $K_{a}$ do

9: if $q_{a}^{k}>0$ then

10: The aggregator computes
$\tilde{R}_{a}^{pf}(q_{a}^{k-})=R_{a}^{v}(q_{a}^{k-})-C_{a}(q_{a}^{k-})$,
where $R_{a}^{v}(q_{a}^{k-})$ is the left-handed limit of the revenue at
$q_{a}^{k}$ computed by (20) and $C_{a}(q_{a}^{k-})$ is the left-handed limit
of the cost at $q_{a}^{k}$ computed by solving Problem CO as in Proposition
3(b);

11: end if

12: The aggregator computes
$\tilde{R}_{a}^{pf}(q_{a}^{k+})=R_{a}^{v}(q_{a}^{k+})-C_{a}(q_{a}^{k+})$,
where $R_{a}^{v}(q_{a}^{k+})$ is the right-handed limit of the revenue at
$q_{a}^{k}$ computed by (20) and $C_{a}(q_{a}^{k+})$ is the right-handed limit
of the cost at $q$ computed by solving Problem CO as in Proposition 3(b);

13: if $\tilde{R}_{a}^{pf}(q_{a}^{k-})>0$ then

14: The aggregator executes Subroutine 3.1 and computes
$\big{(}\hat{q}^{l},\varepsilon^{j},P(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j})\big{)}=QZ^{int}(q_{a}^{k});$

15: break;

16: end if

17: if $\tilde{R}_{a}^{pf}(q_{a}^{k-})=0$ then

18: The aggregator executes Subroutine 3.2 and computes
$\big{(}\hat{q}^{l},\varepsilon^{j},P(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j})\big{)}=QZ^{-}(q_{a}^{k});$

19: break;

20: end if

21: if $\tilde{R}_{a}^{pf}(q_{a}^{k+})\geq 0$ then

22: The aggregator executes Subroutine 3.3 and computes
$\big{(}\hat{q}^{l},\varepsilon^{j},P(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j})\big{)}=QZ^{+}(q_{a}^{k});$

23: break;

24: end if

25: end for

26: output: ($\hat{q}^{l},\varepsilon^{j},P(\hat{q}^{l},\varepsilon^{j}),\
X(\hat{q}^{l},\varepsilon^{j})$);

Subroutine 3.1: $QZ^{int}(q_{a}^{k})$ for Case 1

1: input: price $q_{a}^{k}$;

2: initialization: set iteration index $j=0$, $\varepsilon^{0}>0$, absolute
error $err_{3}>0$ and $err_{4}>0$;

3: The aggregator computes the near-LNP price $\hat{q}^{l}$ (which is also the
LNP price ${q}^{l}$ in the limiting case) between the adjacent threshold
prices $q_{a}^{k-1}$ and $q_{a}^{k}$ as follows:
$\displaystyle\hat{q}^{l}=\frac{q_{a}^{k-1}\tilde{R}_{a}^{pf}(q_{a}^{k-})-q_{a}^{k}\tilde{R}_{a}^{pf}(q_{a}^{(k-1)+})}{\tilde{R}_{a}^{pf}(q_{a}^{k-})-\tilde{R}_{a}^{pf}(q_{a}^{(k-1)+})};$

4: repeat

5: $j=j+1$;

6: $\varepsilon^{j}=\varepsilon^{j-1}/10$;

7: The aggregator executes Subroutine 1 and computes
$\left(R_{a}^{pf}(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j}),P(\hat{q}^{l},\varepsilon^{j})\right)={CU}(\hat{q}^{l},\varepsilon^{j});$

8: until $\mid{R_{a}^{pf}(\hat{q}^{l},\varepsilon^{j})\mid}\leq err_{4};$

9: output:
($\hat{q}^{l},\varepsilon^{j},P(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j})$);

Subroutine 3.2: $QZ^{-}(q_{a}^{k})$ for Case 2

1: input: price $q_{a}^{k}$;

2: initialization: set iteration index $j=0$, $\varepsilon^{0}>0$, absolute
error $err_{3}>$ and $err_{4}>$;

3: The aggregator calculates the slope $slp(q_{a}^{k})$ of $R_{a}^{pf}(q)$
over the threshold price interval $(q_{a}^{k-1},q_{a}^{k})$:
$slp(q_{a}^{k})=\frac{\tilde{R}_{a}^{pf}(q_{a}^{k-})-\tilde{R}_{a}^{pf}(q_{a}^{(k-1)+})}{q_{a}^{k}-q_{a}^{k-1}};$

4: The aggregator computes the near-LNP price $\hat{q}^{l}$ (smaller than
$q_{a}^{k}$) such that $R_{a}^{pf}(\hat{q}^{l})$ is close enough to zero
(within the absolute error $err_{3}$):
$\displaystyle\hat{q}^{l}=q_{a}^{k}-\frac{err_{3}}{slp(q_{a}^{k})};$

5: repeat

6: $j=j+1$;

7: $\varepsilon^{j}=\varepsilon^{j-1}/10$;

8: The aggregator executes Subroutine 1 and computes
$\left(R_{a}^{pf}(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j}),P(\hat{q}^{l},\varepsilon^{j})\right)={CU}(\hat{q}^{l},\varepsilon^{j});$

9: until $\mid
R_{a}^{pf}(\hat{q}^{l},\varepsilon^{j})-R_{a}^{pf}(\hat{q}^{l})\mid\leq
err_{4};$

10: output:
($\hat{q}^{l},\varepsilon^{j},P(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j})$);

Subroutine 3.3: $QZ^{+}(q_{a}^{k})$ for Case 3

1: input: price $q_{a}^{k}$;

2: initialization: set iteration index $j=0$, $\varepsilon^{0}>0$, absolute
error $err_{3}>0$ and $err_{4}>0$;

3: The aggregator calculates the slope $slp(q_{a}^{k+1})$ of $R_{a}^{pf}(q)$
over the threshold price interval $(q_{a}^{k},q_{a}^{k+1})$:
$slp(q_{a}^{k+1})=\frac{\tilde{R}_{a}^{pf}(q_{a}^{(k+1)-})-\tilde{R}_{a}^{pf}(q_{a}^{k+})}{q_{a}^{k+1}-q_{a}^{k}};$

4: The aggregator computes the near-LNP price $\hat{q}^{l}$ (larger than
$q_{a}^{k}$) such that $R_{a}^{pf}(\hat{q}^{l})$ is close enough to
$\tilde{R}_{a}^{pf}(q_{a}^{k+})$ (within the absolute error $err_{3}$):
$\displaystyle\hat{q}^{l}=q_{a}^{k}+\frac{err_{3}}{slp(q_{a}^{k+1})};$

5: repeat

6: $j:=j+1$;

7: $\varepsilon^{j}=\varepsilon^{j-1}/10$;

8: The aggregator executes Subroutine 1 and computes
$\left(R_{a}^{pf}(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j}),P(\hat{q}^{l},\varepsilon^{j})\right)={CU}(\hat{q}^{l},\varepsilon^{j});$

9: until $\mid
R_{a}^{pf}(\hat{q}^{l},\varepsilon^{j})-R_{a}^{pf}(\hat{q}^{l})\mid\leq
err_{4};$

10: output:
($\hat{q}^{l},\varepsilon^{j},P(\hat{q}^{l},\varepsilon^{j}),X(\hat{q}^{l},\varepsilon^{j})$);

## Appendix K: Simulation

### A. Benchmark

For the benchmark system for comparison, we consider a system (without virtual
storage) where each user invests in a physical storage product (e.g., Tesla
Powerwall) by himself. He will optimize the (fixed) capacity of the physical
storage, and can only use his own storage during the lifespan of the storage.
Apart from the capacity cost, each user also bears the power rating cost and
the operational cost by himself. User $i$ solves an optimization problem
$\textbf{BM}_{i}$ as follows to minimize his cost and determines the optimal
storage size over the investment phase.

$\displaystyle\min\ $
$\displaystyle\sum_{\omega}\rho^{\omega}\big{[}C_{i}^{e}(\bm{p}_{i}^{\omega,r,u},\bm{p}_{i}^{\omega,ch},\bm{p}_{i}^{\omega,dis})-R_{i}^{r}(\bm{p}_{i}^{\omega,r,u})$
$\displaystyle+{c^{s}\sum_{t\in\mathcal{T}}({p_{i}^{\omega,ch}[t]+p_{i}^{\omega,dis}[t]})\big{]}+\kappa
c_{u}^{x}x_{i}+\kappa c_{u}^{p}p_{i}}$ s.t.
$\displaystyle\text{(3)},\text{(5)},\text{(6)},\text{(7)},~{}\forall\omega\in\Omega;$
$\displaystyle\gamma^{\text{min}}x_{i}\leq
e_{i}^{\omega}[t]\leq{\gamma^{\text{max}}x_{i}},~{}\forall
t\in\mathcal{T^{\prime}},\forall\omega\in\Omega;$ (54) $\displaystyle{0\leq
p_{i}^{\omega,ch}[t],p_{i}^{\omega,dis}[t]\leq p_{i}},\forall
t\in\mathcal{T},~{}\forall\omega\in\Omega;$ (55) var: $\displaystyle
x_{i},p_{i},\\{\bm{p}_{i}^{\omega,r,u},\bm{p}_{i}^{\omega,ch},\bm{p}_{i}^{\omega,dis},\bm{e}_{i}^{\omega},\forall\omega\in\Omega\\},$

where $c_{u}^{x}$ is the unit capacity cost and $c_{u}^{p}$ is the unit power
rating cost. For Tesla Powerwall, customers need to pay about $7000 (including
installation fee) for the storage capacity of 14KWh and an extra inverter
cost[3]. 212121At present, many companies only provide several capacity
choices for consumers. For example, Tesla Powerwall only offers fixed capacity
of 13.5KWh to consumers. In benchmark $\textbf{BM}_{i}$, we let users flexibly
determine the capacity. Compared with user’s Problem
$\textbf{UP}_{i}^{\omega}$ in our virtualization model, in Problem
$\textbf{BM}_{i}$ each user’s decision of storage capacity $x_{i}$ and power
rating $p_{i}$ are fixed for all scenarios as shown in (54) and (55). Apart
from the capacity cost $c_{u}^{x}x_{i}$, each user $i$ also bears the
operational cost incurred by the charge decision $\bm{p}_{i}^{\omega,ch}$ and
discharge decision $\bm{p}_{i}^{\omega,dis}$ as well as the inverter/converter
cost for the power rating $p_{i}$.

### B. Simulation parameter

We consider the lithium-ion battery as the energy storage technology. We use
realistic load data from PG&E Corporation in 2012 [40] to simulate users’ load
profiles, and we use wind speed and solar radiation data from Hong Kong
Observatory[41] to calculate users’ renewable generations. To illustrate our
storage virtualization model, we simulate a system with three users of
different types, and choose 7 typical scenarios to approximate the original
scenario set.

We use the electricity price data from [28], which is a monthly demand charge
tariff. However, note that in users’ models, we consider a daily billing
period and thus adopt the daily demand charge tariff. Under a monthly demand
charge tariff, the peak power consumption is charged only once over 28-31 days
(depending on which month is concerned). However, in a daily demand charge
tariff, the peak power consumption is charged on each day, and thus 28-31
times in a month. Therefore, we need to scale down the monthly peak price to
the daily peak price, so that users pay the same or similar bill under the
daily and monthly demand charge tariff in each month. We note that it is
nontrivial regarding how to scale down and obtain an accurate daily peak price
for users since different users have different daily and monthly peak load
patterns. Instead, for an illustrative purpose, we calculate the peak price of
the daily demand charge tariff in our simulations as follows. Specifically,

Figure 12: Type-3 user’s load in a month.

we denote the peak price for the monthly demand charge tariff and daily demand
charge tariff by $\pi_{p}^{M}$ and $\pi_{p}^{D}$, respectively. We denote the
peak power in one month of $M$ days (e.g., 30 days) by $p^{m}$ and the peak
power in each day $d$ by $p^{d}$, where $d=1,2,3,...,M$. In order to scale
$\pi_{p}^{M}$ to $\pi_{p}^{D}$, we argue that users maintain the same total
electricity bill under the daily and monthly demand charge, i.e.,

$\sum_{d=1}^{M}p^{d}\pi_{p}^{D}=p^{m}\pi_{p}^{M}.$

Thus, we obtain the scaled-down peak price of the daily demand charge tariff
as

$\pi_{p}^{D}=\pi_{p}^{M}\frac{p^{m}}{\sum_{d=1}^{M}p^{d}}.$

In practice, the daily peak power $p^{d}$ of users in one month can be similar
across days.222222In Figure 12, we show the hourly load of Type-3 user in the
month of 2012 January (with total 744 hours), where the peak power of each day
is highlighted in red circle. We can see that the daily peak power is similar,
and is around 1.75KW. Therefore, for simplicity and for the illustrative
purpose, we consider the case where each day’s peak power is the same, so we
let $p^{d}=p^{m}$ for all $d$. By choosing $M=30$, we then obtain the scaling
factor $\frac{p^{m}}{\sum_{d=1}^{M}p^{d}}=1/30$. Thus, we obtain the peak
price as $\pi_{p}^{D}=0.4$$/KW for the _daily_ demand charge tariff based on
the realistic peak price $\pi_{p}^{M}=10.8$$/KW of the _monthly_ demand charge
tariff. Furthermore, although we choose a daily operational horizon in our
model, our framework can be extended to different timescales of the
operational horizon and billing cycles, e.g., one month.

Therefore, we scale down and obtain the daily peak price $\pi_{p}$ at
$0.4$$/kW. The energy charge price $\pi_{b}$ is $0.03\$$/kWh [28]. We set the
renewable energy selling price $\pi_{s}$ at $0.01\$$/kWh. For the storage
cost, we set the capacity cost $c^{X}$=160$/KWh, power rating cost
$c^{P}$=55$/KW, 232323The lithium-ion battery production cost is estimated
from 145$/KWh to 227$/KWh [46][47][48] in 2017.  and the operational cost
$c^{s}=0.001$$/kWh[35]. For other technical parameters of the storage, we set
efficiency rates $\eta^{c}=\eta^{d}=\eta_{a}^{c}=\eta_{a}^{d}=0.95$, and the
effective capacity rate $\gamma^{\text{min}}=0.9,\gamma^{\text{max}}=1$[2].
For the additional resources, we set $c_{a}^{ch}=0$ and $c_{a}^{dis}=0.1$$/KWh
[49]. We choose a sufficiently small penalty coefficient $\varepsilon=3\times
10^{-7}\$/(\text{kWh})^{2}$.

In addition, our optimization framework is applicable to any parameter
configuration of the system. Interested users can also use other parameter
values of the system to test the performance of our framework.

## Appendix L: The daily capital recovery factor $\kappa$

The aggregator’s investment cost is

$\displaystyle C_{a}^{cap}(X,P)=\kappa c^{X}X+\kappa c^{P}P.$

Note that in practice $c^{X}$ and $c^{P}$ are given as unit costs over an
investment phase of several years. Here, we scale them into the operational
horizon of one day by the daily capital recovery factor $\kappa$. For the
factor $\kappa$, we first calculate the present value of an annuity (a series
of equal annual cash flows) with the annual interest rate $r$, and then we
divide the annuity equally to each day. This leads to the formulation of the
factor $\kappa$ as follows[50],

$\displaystyle\kappa=\frac{r(1+r)^{y}}{(1+r)^{y}-1}\cdot\frac{1}{Y_{d}},$ (56)

where $y$ is the number of years over the total time horizon, and $Y_{d}$ is
the number of days (e.g., 365) in one year.

## Appendix M: Simulation of the case without wind energy

We conduct the simulations where one user has no renewable generation and the
other two users have solar energy. It can be another common case in practice
especially for some areas that are unsuitable to install wind turbines. Recall
that we consider three types of users with seven scenarios in the simulation
part of Section V: Type-1 user has the wind energy while Type-2 and Type-3
users have solar energy. In this part, we conduct the simulation for the case
where Type-1 user has no wind energy while Type-2 and Type-3 users have solar
energy. For the parameter setup, we keep the parameters of the load and
renewable generation the same as the setup in Section V except that Type-1
user has no wind energy. In Figure 13, we show the seven scenarios of load and
renewable generation for Type-1 user (without wind), Type-2 user and Type-3
user, respectively.

We can show that at the optimal-profit price $q^{\star}$, the invested
physical capacity is reduced by 35% compared with the sold virtual capacity.
In Figures 14 and 14, we also show users’ cost reduction in our model compared
with the benchmarks (where users install their own physical storage at the
production cost $c^{p}$ and the market price $c^{r}$, respectively) under the
price $q^{l}$ and $q^{\star}$. Type-2 and Type-3 users’ costs can be reduced
by up to 20.4%, while Type-1 user’s cost can only be reduced by up to 3.1%.
This shows that even though there is no wind energy considered for Type-1
user, our model can still work well for significantly reducing the invested
physical storage capacity and reducing the cost of Type-2 and Type-3 users
with solar energy. Note that in Section V where Type-1 user has the wind
generation, Type-1 user’s cost can be reduced by up to 34.7% (instead of 3.1%
when he has no wind generation). This shows that users without renewable
energy may benefit less compared to users with renewable energy in our
framework.

Figure 13: 7 typical load and renewable generation scenarios for (a) Type-1
user (without wind); (b) Type-2 user; (c) Type-3 user.

Figure 14: (a) Cost reduction at $q^{l}$; (b) Cost reduction at $q^{\star}$.

## Appendix N: Discussion of the uncertainty of the load and renewable
generation

Next we discuss the impact of uncertainties of renewable generations and loads
on the virtual storage operation. In our model, users make purchase decisions
on the virtual capacity as well as the charge and discharge decisions at the
beginning of each day, based on the prediction of their loads and renewable
generations for the next day. There have been extensive studies on the
prediction of wind power generation [8], solar power generation [9], and power
consumption [51]. Since the focus of our work is on the design of the virtual
storage sharing framework, we have initially chosen to assume that users can
perfectly predict their renewable generations and loads. Here we further
provide the discussions about the impact of uncertainties in the day-ahead
prediction on the users’ purchase of virtual capacities and their virtual
storage operation.

To study the impact of uncertainties in loads and renewable generations, we
incorporate prediction errors in the simulations for users’ load and renewable
profiles. Recall that we consider three types of users with seven scenarios in
the simulation part of Section V. For the impact of uncertainties in loads and
renewable profiles, we focus on one type of users and consider one scenario of
the load and renewable profiles as the benchmark of the day-ahead prediction.
Specifically, we choose Type-2 user of Scenario 1, whose (day-ahead predicted)
load and renewable profiles are shown as the solid lines in Figure 15. We
assume that the prediction error $err^{l}$ of the load profile is
characterized by the uniform distribution $U(-\beta P_{i}^{\omega,l}[t],\beta
P_{i}^{\omega,l}[t])$ over intervals $[-\beta P_{i}^{\omega,l}[t],\beta
P_{i}^{\omega,l}[t]]$, where the coefficient $\beta$ denotes the maximum
deviation range of the error $err^{l}$. Thus, in the real-time operation, the
load profile is characterized by $P_{i}^{\omega,l}[t]+err^{l}$ when
considering the prediction error. Similarly, for the renewable generation, we
assume that the prediction error $err^{r}$ is characterized by the uniform
distribution $U(-\beta P_{i}^{\omega,r}[t],\beta P_{i}^{\omega,r}[t])$. Then,
in the real-time operation, the renewable generation is characterized by
$P_{i}^{\omega,r}[t]+err^{r}$ when considering the prediction error. For the
real-time load and renewable profiles, we randomly generate 50 realizations of
the random load $P_{i}^{\omega,l}[t]+err^{l}$ and renewable generation
$P_{i}^{\omega,r}[t]+err^{r}$ with $\beta=10\%$, respectively, the maximum
deviation ranges (away from day-ahead prediction) of which are depicted by the
dashed lines in Figure 15.

We first solve the users’ optimization problem $\textbf{UPP}_{i}^{\omega}$
with the day-ahead prediction load and renewable generation, and obtain the
optimal decisions of the charge/discharge and the purchased virtual capacities
as a benchmark. Then, given 50 realizations of the real-time load and
renewable profiles, we also solve the users’ optimization problem
$\textbf{UPP}_{i}^{\omega}$, and obtain the corresponding optimal decisions of
the charge/discharge and the purchased capacity. We regard the optimal
decisions under the day-ahead predicted profiles as the benchmark. For the
optimal decisions under 50 realizations of the real-time load and renewable
profiles, we show their maximum deviation from the benchmark in Figure 16.

Specifically, for users’ optimal decisions, we assume that the aggregator
chooses the optimal-profit price $p^{\star}$ computed under the day-ahead
prediction. In Figure 16, we show the optimal charge and discharge decision at
the price $p^{\star}$ under the day-ahead prediction as the black curve (where
the positive values represent charge and the negative values represent
discharge), and we show the maximum deviations of the optimal decisions under
the prediction error away from the decisions of the day-ahead prediction as
the blue error bars. In Figure 16, we show the optimal purchased virtual
capacities from the price zero to price $p^{\star}$ under the day-ahead
prediction as the black curve, and we show the maximum deviations under the
prediction error away from the day-ahead prediction as the blue error bars.

In Figure 16, we can see that, when considering the prediction error, the
whole-day optimal charge and discharge decisions will not deviate more than
0.3KW compared with the optimal decisions under the day-ahead prediction. As
shown in Figure 16, the optimal purchased capacity considering the prediction
error will not deviate more than $9\%$ compared with the decisions under the
day-ahead prediction. Furthermore, the cost of Type-2 user under the
prediction error of $\beta=10\%$ is within the range $[47.5,~{}53.1]$ cents,
with a maximum deviation of 6% from the minimized cost 49.87 cents under the
day-ahead prediction.

Figure 15: Type-2 user’s load and renewable profiles.

Figure 16: (a) Type-2 user’s optimal charge and discharge decisions; (b)
Type-2 user’s optimal purchased virtual capacity.

## Appendix O: Diversity of users’ load profiles in a localized region

We refer to the residential load data from [52] and demonstrate that users’
load profiles can be significantly diverse in a localized region.

Figure 17: Load profiles of four households the Mueller community of
Austin,Texas U.S. [10].

Specifically, we acquired residential load data from [52], which is provided
by the Pecan Street Smart Grid Demonstration Program [10]. Over 1000
households in the Mueller community of Austin, Texas, U.S., participated in
the program and shared their electricity consumption data with the project. We
take the hourly load profiles of four different households on January 1, 2017
as an illustrative example to show diverse electricity consumption in Figure
17. We see that Household 1’s load profile is flat across the whole day, but
Household 2 has a peak load at late night (around 22:00). In contrast, the
peak load of Household 3 appears in the evening (around 17:00), and the peak
load of Household 4 occurs around noon (around 11:00).

Thus, we see that the electricity consumption of different households in one
community can be considerably different. Such different load profiles lead to
diverse charge and discharge decisions of users, which can potentially reduce
the required capacity of the physical storage and save the investment cost as
demonstrated in the simulations in Section V of the main text of our paper.

Figure 18: Net load of (a) Type-1 user; (b) Type-2 user; (c) Type-3 user; (d)
the system.

## Appendix P: Peak reduction in the system

Even though we do not explicitly consider minimizing the system peak load as
the objective in our model, we numerically demonstrate that the virtual charge
and discharge of users not only helps users cut their electricity bill but
also leads to a peak load reduction in the system

For the purpose of illustration, we depict the results of the peak-load
reduction of individual users and the peak-load reduction of the system in
Figure 18 of this response letter, considering a three-user system. In Figure
18, 18 and 18, we show three types of users’ original net daily load (i.e.,
load minus the renewable generation) compared with the net load after
utilizing the virtual storage (where the blue dashed curve corresponds to the
optimal-profit (OP) price $q^{\star}$ and the red dotted curve corresponds to
the lowest-nonnegative-profit (LNP) price $q^{l}$). We see that each user’s
peak net-load is shaved significantly by more than 40%.242424Note that for
Type-2 and Type-3 users, their decisions of virtual storage are the same at OP
price $q^{\star}$ and at LNP price $q^{l}$, and thus the net load curves
overlap with each other. In Figure 18, we show the aggregated original system
net load (i.e., aggregated load minus the aggregated renewable generation of
all three users) compared with the aggregated system net load after utilizing
the virtual storage (where the blue dashed curve corresponds to the OP price
$q^{\star}$ and the red dotted curve corresponds to the LNP price $q^{l}$). We
see that the aggregated peak load of the system in one day is reduced by 41.7%
at $q^{\star}$ and 44.0% at $q^{l}$, compared with the system peak load
without virtual storage. Thus, the results shown in Figure 18 demonstrate that
our model not only reduces users’ peak load, but also reduces the system peak,
and therefore benefits the whole system.

## References

  * [1] D. Zhao, H. Wang, J. Huang, and X. Lin, “Pricing-based energy storage sharing and virtual capacity allocation,” in _Communications (ICC), 2017 IEEE International Conference on_ , May 2017, pp. 1–6.
  * [2] X. Luo, J. Wang, M. Dooner, and J. Clarke, “Overview of current development in electrical energy storage technologies and the application potential in power system operation,” _Applied Energy_ , vol. 137, pp. 511–536, 2015.
  * [3] “Tesla powerwall.” [Online]. Available: http://www.energysage.com/solar/solar-energy-storage/tesla-powerwall-home-battery
  * [4] Q. Zhang, L. Cheng, and R. Boutaba, “Cloud computing: state-of-the-art and research challenges,” _Journal of internet services and applications_ , vol. 1, no. 1, pp. 7–18, 2010.
  * [5] E. Dickinson, _Solar Energy Technology Handbook: 0_. CRC Press, 2018.
  * [6] A. C. Orrell, N. A. Foster, J. S. Homer, and S. L. Morris, “2015 distributed wind market report,” Pacific Northwest National Lab.(PNNL), Richland, WA (United States), Tech. Rep., 2016.
  * [7] H. Wang and J. Huang, “Joint investment and operation of microgrid,” _IEEE Trans. on Smart Grid_ , vol. 8, no. 2, pp. 833–845, March 2017.
  * [8] A. M. Foley, P. G. Leahy, A. Marvuglia, and E. J. McKeogh, “Current methods and advances in forecasting of wind power generation,” _Renewable Energy_ , vol. 37, no. 1, pp. 1–8, 2012.
  * [9] R. H. Inman, H. T. Pedro, and C. F. Coimbra, “Solar forecasting methods for renewable energy integration,” _Progress in energy and combustion science_ , vol. 39, no. 6, pp. 535–576, 2013.
  * [10] “Pecan street grid demonstration program. final technology performance report,” Pecan Street Project, Inc., Austin, TX (United States), Tech. Rep., 2015\.
  * [11] D. Zhao, H. Wang, J. Huang, and X. Lin, “Technical report for ”virtual energy storage sharing and capacity allocation” (2019).” [Online]. Available: https://arxiv.org/abs/1907.02005
  * [12] J. D. Craig, “Motivations for market restructuring: evidence from us electricity deregulation,” _Energy Economics_ , vol. 60, pp. 162–167, 2016\.
  * [13] E. Parliament, “Understanding electricity markets in the eu,” 2016.
  * [14] J. Jin and Y. Xu, “Optimal storage operation under demand charge,” _IEEE Trans. on Power Systems_ , vol. 32, no. 1, pp. 795–808, Jan 2017.
  * [15] C. O. Adika and L. Wang, “Non-cooperative decentralized charging of homogeneous households’ batteries in a smart grid,” _IEEE Trans. on Smart Grid_ , vol. 5, no. 4, pp. 1855–1863, July 2014.
  * [16] I. Atzeni, L. G. Ordóñez, G. Scutari, D. P. Palomar, and J. R. Fonollosa, “Demand-side management via distributed energy generation and storage optimization,” _IEEE Trans. on Smart Grid_ , vol. 4, no. 2, pp. 866–876, June 2013.
  * [17] Z. Wang, C. Gu, F. Li, P. Bale, and H. Sun, “Active demand response using shared energy storage for household energy management,” _IEEE Trans. on Smart Grid_ , vol. 4, no. 4, pp. 1888–1897, Dec 2013.
  * [18] J. Yao and P. Venkitasubramaniam, “Stochastic games of end-user energy storage sharing,” in _Proc.of IEEE CDC_ , Dec 2016, pp. 4965–4972.
  * [19] ——, “Optimal end user energy storage sharing in demand response,” in _IEEE SmartGridComm_ , Nov 2015, pp. 175–180.
  * [20] C. P. Mediwaththe, E. R. Stephens, D. B. Smith, and A. Mahanti, “A dynamic game for electricity load management in neighborhood area networks,” _IEEE Trans. on Smart Grid_ , vol. 7, no. 3, pp. 1329–1336, May 2016.
  * [21] W. Tushar, B. Chai, C. Yuen, S. Huang, D. B. Smith, H. V. Poor, and Z. Yang, “Energy storage sharing in smart grid: A modified auction-based approach,” _IEEE Transactions on Smart Grid_ , vol. 7, no. 3, pp. 1462–1475, May 2016\.
  * [22] H. Chen, Y. Yu, Z. Hu, H. Luo, C. Tan, and R. Rajagopal, “Energy storage sharing strategy in distribution networks using bi-level optimization approach,” in _2017 IEEE Power Energy Society General Meeting_ , July 2017, pp. 1–5.
  * [23] J. Liu, N. Zhang, C. Kang, D. S. Kirschen, and Q. Xia, “Decision-making models for the participants in cloud energy storage,” _IEEE Trans. on Smart Grid_ , vol. PP, no. 99, pp. 1–1, 2017.
  * [24] W. Wei, F. Liu, and S. Mei, “Energy pricing and dispatch for smart grid retailers under demand response and market price uncertainty,” _IEEE Transactions on Smart Grid_ , vol. 6, no. 3, pp. 1364–1374, May 2015.
  * [25] H. Wang and J. Huang, “Incentivizing energy trading for interconnected microgrids,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 4, pp. 2647–2657, July 2018.
  * [26] F. Lezama, J. Soares, P. Hernandez-Leal, M. Kaisers, T. Pinto, and Z. M. A. do Vale, “Local energy markets: Paving the path towards fully transactive energy systems,” _IEEE Transactions on Power Systems_ , pp. 1–1, 2018.
  * [27] L. Wood, R. Hemphill, J. Howat, R. Cavanagh, S. Borenstein, J. Deason, and L. Schwartz, “Recovery of utility fixed costs: Utility, consumer, environmental and economist perspectives,” Lawrence Berkeley National Lab, Berkeley, CA, Tech. Rep., 2016.
  * [28] R. Hledik, “Rediscovering residential demand charges,” _The Electricity Journal_ , vol. 27, no. 7, pp. 82–96, 2014.
  * [29] S. Zhao, X. Lin, and M. Chen, “Peak-minimizing online ev charging: Price-of-uncertainty and algorithm robustification,” in _IEEE INFOCOM_ , 2015, pp. 2335–2343.
  * [30] R. L. Fares and M. E. Webber, “The impacts of storing solar energy in the home to reduce reliance on the utility,” _Nature Energy_ , vol. 2, no. 2, p. 17001, 2017.
  * [31] M. Cheng, S. S. Sami, and J. Wu, “Virtual energy storage system for smart grids,” _Energy Procedia_ , vol. 88, pp. 436–442, 2016.
  * [32] ——, “Benefits of using virtual energy storage system for power system frequency response,” _Applied energy_ , vol. 194, pp. 376–385, 2017.
  * [33] J. Kwac, J. Flora, and R. Rajagopal, “Lifestyle segmentation based on energy consumption data,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 4, pp. 2409–2418, July 2018.
  * [34] J. A. Jardini, C. M. V. Tahan, M. R. Gouvea, S. U. Ahn, and F. M. Figueiredo, “Daily load profiles for residential, commercial and industrial low voltage consumers,” _IEEE Transactions on Power Delivery_ , vol. 15, no. 1, pp. 375–380, Jan 2000.
  * [35] H. Wang and J. Huang, “Incentivizing energy trading for interconnected microgrids,” _IEEE Trans. on Smart Grid_ , vol. PP, no. 99, pp. 1–13, 2016\.
  * [36] Y. Guo and Y. Fang, “Electricity cost saving strategy in data centers by using energy storage,” _IEEE Trans. on Parallel and Distributed Systems_ , vol. 24, no. 6, pp. 1149–1160, 2013.
  * [37] E. A. Ok, _Real analysis with economic applications_. Princeton University Press, 2007.
  * [38] G. Dantzig, _Linear programming and extensions_. Princeton university press, 2016.
  * [39] A. Altman and J. Gondzio, “Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization,” _Optimization Methods and Software_ , vol. 11, no. 1-4, pp. 275–302, 1999\.
  * [40] “Understanding energy use and prices,.” [Online]. Available: https://www.pge.com/tariffs/energy_use_prices.shtml
  * [41] “Summary of wind and solar meteorological data in hong kong.” [Online]. Available: http://jianwei.ie.cuhk.edu.hk/energy.html
  * [42] U. D. of Energy, “The rise of community wind power in japan: Enhanced acceptance through social innovation,” 2012.
  * [43] Y. Maruyama, M. Nishikido, and T. Iida, “Community wind benefits,” _Energy Policy_ , vol. 35, no. 5, pp. 2761–2769, 2007.
  * [44] “Data of vritual storage sharing.” [Online]. Available: https://drive.google.com/open?id=1i5FCmfiwU33Sv0iogQWGASik0Bykzfjl
  * [45] A. B. Berkelaar, K. Roos, and T. Terlaky, “The optimal set and optimal partition approach to linear and quadratic programming,” in _Advances in Sensitivity Analysis and Parametic Programming_. Springer, 1997, pp. 159–202.
  * [46] “Chevy bolt production confirmed for 2016.” [Online]. Available: http://www.hybridcars.com/chevy-bolt-production-confirmed-for-2016
  * [47] “Electric vehicle battery cost dropped 80% in 6 years down to $227/kwh.” [Online]. Available: https://electrek.co/2017/01/30/electric-vehicle-battery-cost-dropped-80-6-years-227kwh-tesla-190kwh
  * [48] “Lithium-ion battery costs and markets.” [Online]. Available: https://data.bloomberglp.com/bnef/sites/14/2017/07/BNEF-Lithium-ion-battery-costs-and-market.pdf
  * [49] P. S. E. Inc., “Costs of utility distributed generators, 1-10 mw: Twenty-four case studies,” EPRI, Palo Alto, CA, and Cooperative Research Network, Arlington, VA: 2003. 1007760.
  * [50] H. Pandžić, Y. Wang, T. Qiu, Y. Dvorkin, and D. S. Kirschen, “Near-optimal method for siting and sizing of distributed storage in a transmission network,” _IEEE Trans. on Power Systems_ , vol. 30, no. 5, pp. 2288–2300, Sept 2015.
  * [51] D. Lachut, N. Banerjee, and S. Rollins, “Predictability of energy use in homes,” in _Green Computing Conference (IGCC), 2014 International_. IEEE, 2014, pp. 1–10.
  * [52] “Dataport from pecan street.” [Online]. Available: https://dataport.cloud.

|  Dongwei Zhao Dongwei Zhao (S’16) received the B.E. degree from the College
of Electrical Engineering, Zhejiang University, Hangzhou, China, in 2015. He
is currently pursuing the Ph.D. degree with the Department of Information
Engineering, The Chinese University of Hong Kong. His research areas include
optimization and game theory of power and energy systems.
---|---
|  Hao Wang Hao Wang (M’16) is a Postdoctoral Scholar at Stanford University.
He received his Ph.D. from The Chinese University of Hong Kong and has been a
Washington Research Foundation Innovation Fellow at University of Washington,
Seattle. His main research interests are in the optimization, machine
learning, and data analytics of power and energy systems. More information at
https://web.stanford.edu/~haowang6/.
---|---
|  Jianwei Huang Jianwei Huang (F’16) is a Presidential Chair Professor and
the Associate Dean of the School of Science and Engineering, The Chinese
University of Hong Kong, Shenzhen. He is also a Professor in the Department of
Information Engineering, The Chinese University of Hong Kong. He is the co-
author of 9 Best Paper Awards, including IEEE Marconi Prize Paper Award in
Wireless Communications 2011. He has co-authored six books, including the
textbook on ”Wireless Network Pricing”. He has served as the Chair of IEEE
ComSoc Cognitive Network Technical Committee and Multimedia Communications
Technical Committee. He has been an IEEE Fellow, an IEEE ComSoc Distinguished
Lecturer, and a Clarivate Analytics Highly Cited Researcher. More information
at http://jianwei.ie.cuhk.edu.hk/.
---|---
|  Xiaojun Lin (S’02 M’05 SM’12 F’17) received his B.S. from Zhongshan
University, Guangzhou, China, in 1994, and his M.S. and Ph.D. degrees from
Purdue University, West Lafayette, IN, in 2000 and 2005, respectively. He is
currently a Professor of Electrical and Computer Engineering at Purdue
University. Dr. Lin’s research interests are in the analysis, control and
optimization of large and complex networked systems, including both
communication networks and power grid. He received the IEEE INFOCOM 2008 best
paper and 2005 best paper of the year award from Journal of Communications and
Networks. He received the NSF CAREER award in 2007. He was the Workshop co-
chair for IEEE GLOBECOM 2007, the Panel co-chair for WICON 2008, the TPC co-
chair for ACM MobiHoc 2009, the Mini-Conference co-chair for IEEE INFOCOM
2012, and General co-chair for ACM e-Energy 2019. He is currently serving as
an Area Editor for (Elsevier) Computer Networks Journal, and has served as an
Associate Editor for IEEE/ACM Transactions on Networking and a Guest Editor
for (Elsevier) Ad Hoc Networks journal.
---|---