# Laser-power consumption of soliton formation in a bidirectional Kerr
resonator

Jizhao Zang1,2∗, Su-Peng Yu1,2, Haixin Liu1,2, Yan Jin1,2, Travis C. Briles1,
David R. Carlson1,3, and Scott B. Papp1,2

1Time and Frequency Division, National Institute of Standards and Technology,
Boulder, CO, USA
2Department of Physics, University of Colorado, Boulder, CO, USA
2Octave Photonics, Louisville, CO, USA

∗E-mail<EMAIL_ADDRESS>
> Laser sources power extreme data transmission as well as computing
> acceleration, access to ultrahigh-speed signaling, and sensing for
> chemicals, distance, and pattern recognition. The ever-growing scale of
> these applications drives innovation in multi-wavelength lasers for
> massively parallel processing. We report a nanophotonic Kerr-resonator
> circuit that consumes the power of an input laser and generates a soliton
> frequency comb at approaching unit efficiency. By coupling forward and
> backward propagation, we realize a bidirectional Kerr resonator that
> supports universal phase matching but also opens excess loss by double-sided
> emission. Therefore, we induce reflection of the resonator’s forward,
> external-coupling port to favor backward propagation, resulting in
> efficient, one-sided soliton formation. Coherent backscattering with
> nanophotonics provides the control to put arbitrary phase-matching and
> efficient laser-power consumption on equal footing in Kerr resonators. In
> the overcoupled-resonator regime, we measure 65% conversion efficiency of a
> 40 mW input pump laser, and the nonlinear circuit consumes 97% of the pump,
> generating the maximum possible comb power. Our work opens up high-
> efficiency soliton formation in integrated photonics, exploring how energy
> flows in nonlinear circuits and enabling laser sources for advanced
> transmission, computing, quantum sensing, and artificial-intelligence
> applications.

##

Nonlinear optics is a ubiquitous phenomenon in science and technology that
converts light from one wavelength to another and enables the synthesis of new
optical fields. Development in nonlinear optics has advanced a host of
applications, ranging from optical signal processing (?) to ultra-short pulsed
lasers (?) and optical parametric amplifiers (?, ?, ?) to quantum-state
generation and manipulation (?, ?, ?, ?, ?). Invariably, nonlinear effects are
relatively weak, and harnessing their benefits stimulates discovery and
advancements in engineering the control of light-matter interactions. Since
the third-order Kerr-effect nonlinearity involving intensity-dependent
refractive index occurs in practically all materials, it is widely applied in
scalable technologies like integrated photonics (?, ?, ?, ?, ?).

Controlling Kerr nonlinearity in microresonators is a central goal for
integrated photonics. The high quality factor and small mode volume of
microresonators make optical parametric oscillation possible with a threshold
power much less than one milliwatt (?), enabling the use with many different
types of lasers. Moreover, adjusting the material and shape of a
microresonator enables phase-matching for four-wave mixing (FWM). There now
exists a fundamental understanding of stability criteria for many types of
nonlinear field states in Kerr microresonators, including parametric
oscillation (?, ?, ?) and soliton microcombs (?, ?, ?, ?). We consider the
latter as self-reinforcing, temporally localized pulses in a resonator that
manifest as a frequency comb at the resonator output. Recently, soliton
stability and generation has been enhanced by coupling forward and backward
propagation in microresonators as an added mechanism to control nonlinear
phase matching. Photonic-crystal resonator (PhCR) microcombs (?, ?, ?)
implement coherent scattering with sub-wavelength nanopatterns to realize a
bidirectional Kerr resonator, and self-injection-locked microcombs exploit
direct laser-resonator feedback (?, ?, ?, ?, ?). These bidirectional systems
can feature universal phase matching in any dispersion regime and spontaneous
soliton formation.

While phase-matching is the essential condition to initiate nonlinear field
states like soliton microcombs, efficiently energizing and harvesting the
soliton comb represents a distinct set of physical optimizations. We define
the pump-to-comb conversion efficiency (CE) as
$\rm{CE}=P_{\rm{comb}}/P_{\rm{pump}}$, where $P_{\rm{pump}}$ is the pump laser
power and $P_{\rm{comb}}$ is the total, emitted comb power. High CE increases
the power per comb mode, which is critical for applications like wavelength-
division-multiplexing (WDM) transceivers (?, ?) and microwave photonic filters
(?). To maximize CE, first, the external coupling rate ($\kappa_{c}$) of the
microresonator must be larger than the internal loss rate ($\kappa_{i}$).
Furthermore, phase-matching constrains the CE of solitons because only limited
ranges of system parameters like laser detuning are allowed. Still, isolated,
high CE soliton states have been reported. Efficient temporal solitons have
been demonstrated in coupled fiber cavities (?) and a microresonator nested in
a fiber loop with gain (?), but those systems are not compatible with chip-
scale integration. Soliton crystals achieve nearly unit CE by exciting several
solitons in one resonator, however, destructive interference reduces the
utility of such ensembles (?). Solitons in anomalous-dispersion resonators
with artificially shifted pump resonance offer a path to high CE (?), but
precise control of the pump laser and auxiliary cavity is required to initiate
the soliton state. Solitons in normal-dispersion resonators (?, ?, ?, ?) can
leverage a long duty cycle for enhanced CE, however soliton phase-matching is
operationally challenging to achieve. Hence, unified control of phase-matching
and CE for high pump-laser consumption in solitons is an outstanding
objective.

Here, we leverage coherent scattering in a nanophotonic Kerr-resonator circuit
to provide universal phase-matching for soliton microcombs and demonstrate
deterministic control of CE up to the unit boundary. We experimentally
generate solitons with a PhCR that features normal dispersion and one ’split
mode’ with coherently coupled forward and backward propagation. Moreover, this
direct phase matching enables a continuous evolution from the empty resonator
state to a backward propagating soliton. However, CE in a bidirectional
resonator cannot exceed 50% because the localized soliton circulates in one
direction while the pump is distributed in both. To consume the forward pump
power and close off excess loss of the double-sided resonator, we place a
reflector on the coupling waveguide. Our nanophotonic soliton circuit emits a
200 GHz train of flattop pulses with up to 65% CE and 97% pump consumption.
Moreover, the nanophotonic circuit enables precise control of the comb
spectrum that we use for alignment to the ITU-T grid with $<2.6$ GHz offset
and compatibility with a distributed feedback (DFB) laser diode on the 193.1
THz channel. Our work demonstrates near unit efficiency comb generation with a
bidirectional Kerr resonator by controlling pump utilization with on-chip
nanophotonics, enabling for example massive scaling in data transmission and
processing.

Figure 1: Soliton formation in the nanophotonic circuit. (A) Schematic of the
nanophotonic Kerr resonator circuit with a PhCR and a waveguide reflector. (B)
Scanning electron microscope images of the nanophotonic structures. (C) The
bidirectional pump in the PhCR supports excitation of solitons in the forward
and backward directions. We tune $\phi$ for constructive interference between
the reflected pump and the backward pump in the PhCR. (D) Power flow in the
circuit analyzed by LLE simulation with $(P_{f}+P_{b})/P_{\rm{pump}}$ (red),
$P_{\rm{comb}}^{f}/P_{\rm{pump}}$ (green), and
$P_{\rm{comb}}^{b}/P_{\rm{pump}}$ (blue). (E) Measured CE versus $\phi$
(circles) and LLE simulation (line). (F) Measured (data trace) and simulated
(line) optical spectrum of a soliton microcomb with a CE of 65%.

Figure 1 introduces the structure of our nanophotonic resonator circuit with
an emphasis on how we achieve unified control of phase-matching and CE through
coherent backscattering. The circuit (Fig. 1A) is composed of: a PhCR in which
the ring width (RW) is modulated with a periodicity of $\pi a/m$, where $a$ is
the PhCR radius and $m$ is the azimuthal mode number; an evanescent coupling
waveguide curved around the PhCR for enhanced escape efficiency; and a
reflector in the coupling waveguide in which the waveguide width is modulated
with a periodicity of $\Lambda=\pi c/\omega_{l}n_{\rm{eff}}$, where c is the
speed of light, $\omega_{l}$ is the laser frequency, and $n_{\rm{eff}}$ is the
effective refractive index. We denote PhCR mode $m$ as the split mode, which
manifests as two resonances split in frequency by the bandgap $\epsilon$ due
to the coherently coupled forward and backward propagation. To energize a
soliton in a normal dispersion PhCR, we tune the pump laser into the lower
frequency resonance of the split mode, which offers optimal phase-matching.
The higher frequency resonance of the split mode does not offer phase-matching
for FWM because the adjacent PhCR modes are all at lower frequencies. Since we
induce nearly complete reflection over a 3 dB bandwidth of 7.6 THz, the
forward-propagating pump laser and backward-propagating soliton microcomb
combine at a circulator and couple to the chip at a common facet. We implement
the circuit on a silicon chip by use of conventional semiconductor processing
with the tantalum pentoxide ($\rm{Ta_{2}O_{5}}$, hereafter tantala) material
platform. We perform nanofabrication on a 75 mm wafer with $<$150 nm feature
resolution in the tantala device layer, creating 50 chips with 100 circuit
variations per chip for experiments (?, ?, ?, ?). Using electron-beam
lithography, we are able to precisely control the circuit dimensions,
especially RW modulation at the nanometer scale. Figure 1B shows scanning
electron microscope images of the nanophotonic circuit elements, emphasizing
the modulation on the inner edge of the PhCR and the width-modulated waveguide
that forms the reflector.

To understand operation of the nanophotonic circuit, we consider the forward
and backward propagating fields; see Fig. 1C. Coherent scattering in the PhCR
creates a bidirectional pump that in principle supports simultaneous
excitation of solitons in the forward and backward directions. However,
backward operation shows higher gain and is favored over forward operation
(?). Hence, the pump laser in the forward direction is not used. The function
of the reflector is to reflect the forward pump field and constructively
interfere with the backward pump field in the PhCR. We characterize the phase
delay $\phi$ of the reflected pump with respect to the PhCR, and
operationally, we control $\phi$ by rotating the PhCR with respect to the
coupling waveguide.

We explore the dynamics and steady-state field solution of the PhCR with the
Lugiato-Lefever equation (LLE) (?), including the forward-backward coupling
and the pump reflector that we implement through coherent scattering (?).
Moreover, the LLE depends on dissipation, group-velocity dispersion (GVD), and
detuning of the pump laser. In particular,

$\begin{split}\frac{\partial E^{f}_{\mu}}{\partial
t}=-(1+i\alpha)E^{f}_{\mu}-i\frac{d_{2}}{2}\mu^{2}E^{f}_{\mu}+i\mathcal{F}(|E^{f}|^{2}E^{f})_{\mu}+2iE^{f}_{\mu}\sum_{j}|E^{b}_{j}|^{2}\\\
+\delta_{\mu 0}(F-\frac{i\epsilon_{\rm{PhC}}}{2}E^{b}_{\mu})\\\ \frac{\partial
E^{b}_{\mu}}{\partial
t}=-(1+i\alpha)E^{b}_{\mu}-i\frac{d_{2}}{2}\mu^{2}E^{b}_{\mu}+i\mathcal{F}(|E^{b}|^{2}E^{b})_{\mu}+2iE^{b}_{\mu}\sum_{j}|E^{f}_{j}|^{2}\\\
+\delta_{\mu
0}(rF-\frac{i\epsilon_{\rm{PhC}}}{2}E^{f}_{\mu})-I_{\Omega}(\mu)\gamma
rE^{f}_{\mu}\end{split}$ (1)

where $E^{f}_{\mu}$ and $E^{b}_{\mu}$ are the forward and backward field
amplitude in mode with index $\mu$ relative to the split mode. All frequency
variables are normalized to the half-linewidth of the resonator $\kappa/4\pi$,
where $\kappa=\kappa_{i}+\kappa_{c}$ is the total loss rate. The variable
$\alpha=2(\omega_{l}-\omega_{0})/\kappa$ is the normalized detuning between
$\omega_{l}$ and the split mode $\omega_{0}$. The parameter $d_{2}$ is the
normalized GVD $D_{2}$ that is derived from the Taylor expansion of the
resonance frequencies around $\omega_{0}$:
$\omega_{\mu}=\omega_{0}+D_{1}\mu+D_{2}\mu^{2}/2+\cdots$, where $D_{1}/2\pi$
is the free-spectral range (FSR). The third and fourth terms of Equation (1)
stand for self-phase modulation and cross-phase modulation, respectively,
where $E^{f}$ and $E^{b}$ are the temporal amplitudes of the intraresonator
field and $\mathcal{F}()_{\mu}$ represents the $\mu$th frequency component of
the Fourier series. Different sources contribute to the pump power in the
forward and backward directions: in the forward direction, the pump power
originates from the external laser on the coupling waveguide and the power
coupled from $E^{b}_{0}$ due to coherent backscattering in the PhCR; in the
backward direction, the pump power originates from the external laser after
travelling through the evanescent coupling region and being reflected,
$E^{f}_{0}$ is outcoupled, reflected and coupled back into the PhCR, and
coherent backscattering of $E^{f}_{0}$ inside the resonator. In the fifth
term, $\delta_{\mu 0}$ is the Kronecker delta ($\delta_{00}=1$, otherwise
$\delta_{\mu 0}=0$), $r=\sqrt{R}\ e^{i\phi}$ is the reflection coefficient of
the waveguide reflector ($R$ is the reflectivity), $F$ is the normalized
amplitude of the pump, and $\epsilon_{\rm{PhC}}=4\pi\epsilon/\kappa$ is the
normalized $\epsilon$. The transfer function from the resonator to the
coupling waveguide is defined as $\gamma=2K/(K+1)$, where $K$ is the coupling
factor. Since the reflector has a finite reflection bandwidth, we apply an
indicator function $I_{\Omega}(\mu)$: if $\mu\in\Omega$, $I_{\Omega}(\mu)=1$,
otherwise $I_{\Omega}(\mu)=0$, where $\Omega$ is a set of $\mu$ that
corresponds to the resonator modes in the reflection band. In LLE simulation,
we use $\Omega=\\{-19,-18,\ldots,18,19\\}$ to model our nanophotonic circuit
and we assume $r$ is the same for $\mu\in\Omega$.

Using the LLE, we simulate power flow in the nanophotonic circuit to
understand CE and pump laser consumption; see Fig. 1D. To quantify the power
flow in both propagating directions on the coupling waveguide, we denote the
comb power and the outcoupled pump power in the forward and backward
directions as $P_{\rm{comb}}^{f}$, $P_{\rm{comb}}^{b}$, $P_{f}$ and $P_{b}$,
respectively. Then, $P_{\rm{comb}}$ naturally becomes
$P_{\rm{comb}}^{f}+P_{\rm{comb}}^{b}$. As a function of $\alpha$, we simulate
$(P_{f}+P_{b})/P_{\rm{pump}}$ (red), $P_{\rm{comb}}^{f}/P_{\rm{pump}}$
(green), and $P_{\rm{comb}}^{b}/P_{\rm{pump}}$ (blue) in the coupling
waveguide. Since the larger gain and the reflector favor backward soliton
formation, we consider $P_{\rm{comb}}^{b}\gg P_{\rm{comb}}^{f}$ and
$P_{\rm{comb}}\approx P_{\rm{comb}}^{b}$ unless otherwise specified. In the
simulation, the optimal reflection phase for maximum CE is $\phi=0$ rad. The
other parameter of this simulation is $K$, which we set to a value of 4.5
based on the likelihood of successfully fabricating a coupler with this value.
As we inspect the simulated behavior of soliton formation in this device in a
similar fashion as we sweep $\alpha$ in the experiment, we observe initial
consumption of the pump laser due to internal losses of the resonator. Then we
observe rapid formation of the soliton, reaching CE and laser consumption of
67% and 97%, respectively. As $\alpha$ exceeds the range over which the
soliton is stable in the resonator, the soliton rapidly de-energizes. These
predicted behaviors are favorable for applications, especially the high CE,
hence we are motivated to explore soliton formation with fabricated versions
of this system.

Our simulations indicate that $\phi$ is the critical parameter to access the
regime of near unit CE in our nanophotonic circuit. Hence, we design a set of
circuits on a common chip in which we systematically and precisely vary
$\phi$; by testing such copies of the circuit, we search for optimum CE. Here,
our test procedure is to scan $\alpha$, record the reflection of the split
mode to characterize $\phi$ dependence of pump interference, and then pump the
lower frequency resonance of the split mode with optimal $\phi$.

Figure 1E presents measurements of CE (circles) in the nanophotonic circuit
alongside the LLE simulation (line) as a function of $\phi$. The agreement of
the simulation and experiment highlights our understanding of the system, our
capability to accurately fabricate specific circuit designs with a set of
$\phi$ values, and robustness of the nonlinear nanophotonic circuit operation.
The variation in $\phi$ results in constructive and destructive interference
between the reflected forward pump by the reflector and backward pump in the
PhCR, leading to different CE in the nanophotonic circuit. Here, we use $K=3$
and the highest CE is 57% in the experimental data, but our simulation
suggests that higher CE is possible with a larger $K$. For example, we observe
a CE of 65% at $K=4.5$ and $P_{\rm{pump}}=40$ mW. In Fig. 1F, we present the
optical frequency spectrum that we obtain with this highest CE device. We tune
$\alpha$ to access this spectrum in which the mode-to-mode power variation is
relatively small. The measured soliton microcomb spectrum (data trace) is in
agreement with the LLE simulation (line). We also observe substantial
consumption of the pump laser in forming solitons, and the remaining pump
power exiting the circuit is only 1 mW.

Figure 2:  Characterization of the nanophotonic circuit structures. (A) Setup
for device characterization. AWG, arbitrary waveform generator; CW laser,
continuous-wave laser; SOA, semiconductor optical amplifier; OSA, optical
spectrum analyzer; PD, photodetector; OSC, oscilloscope. (B) Measured
$\epsilon$ as a function of the modulation amplitude on the inner
($A_{\rm{in}}$, triangles) or outer ($A_{\rm{out}}$, squares) edge of the PhCR
($a=109.5$ $\mu$m, RW = 4 $\mu$m). Inset: transmission of a split mode. (C)
Optimization of the reflection band and $R$: transmission spectra of the
reflector for varied $\Lambda$ and $W_{\rm{PhC}}$. (D) Nanophotonic circuit
reflection as we step $\phi$ in amounts of $0.3\pi$.

We present characterization of the nanophotonic circuit structures that
control CE; see Fig. 2. We use a narrow linewidth laser and a broadband light
source to measure with fine frequency resolution and across the entire
operating bandwidth of the circuit, respectively. We record optical signals
with a standard optical spectrum analyzer and a photodetector with an emphasis
on the reflected power from the circuit; see Fig. 2A.

As we show in Fig. 2B, we implement the azimuthal RW modulation ($A_{\rm{in}}$
and $A_{\rm{out}}$) on either the inner or outer edge of the PhCR, since the
overlap between the interactivity field and the nanopattern permits access to
different ranges of modulation amplitude. With the PhCR, the number of
modulation periods corresponds to twice $m$, the mode number of the PhCR mode
that exhibits coherent backscattering, and the amplitude of the modulation
controls $\epsilon$ of the split mode. We measure $\epsilon$ by a calibrated
laser frequency scan across the split mode, and access to $\epsilon$ from zero
to 5 GHz provides the range that our LLE simulations indicate are needed for
soliton formation. To implement and test the reflector in the coupling
waveguide, we tune $\Lambda$ so that the reflection band is centered on the
split mode, and we optimize the amplitude of the nanopattern for $R>90\%$.
Figure 2C shows our operational procedure to optimize the reflection band as
we vary $\Lambda$ to shift the frequency of maximum reflection. We observe
this shift in the transmission spectra as we step $\Lambda$ in steps
$\Delta\Lambda=10$ nm. This relatively coarse resolution in the nanopattern is
sufficient to align the split mode and peak $R$ of the reflector. We also
characterize the dependence of $R$ on the amplitude of the reflector
nanopattern; see the modulation depth $W_{\rm{PhC}}$ in Fig. 2C. By adjusting
$W_{\rm{PhC}}$ from 750 nm to 850 nm, we vary $R$ from 90% to 97%.

With optimized reflection bandwidth and $R$, we characterize the
intraresonator pump power at the two resonances of the split mode. Figure 2D
shows measurements of reflection from the nanophotonic circuit in which we
step $\phi$ in amounts of $0.3\,\pi$. The data demonstrates interference of
the pump laser, particularly destructive interference of either the higher or
lower frequency component of the split mode. Due to the electric field profile
in the PhCR, the two resonances are out of phase by $\pi$ (?). Operationally,
we vary $\phi$ by rotating the PhCR to change the alignment of the azimuthal
nanopattern and the reflector. This procedure mitigates fabrication
imperfections due to the rotational symmetry of the resonator. Another benefit
of the reflector in our nanophotonic circuit is reduction in thermal
instability from absorbed pump power in the resonator. At the optimal setting
of $\phi$, the higher frequency resonance is not substantially excited and
hence does not contribute to absorption and heating of the PhCR.

Figure 3:  Laser-power consumption and CE in the nanophotonic circuit. (A)
Experimental setup for soliton generation and characterization. EDFA, erbium-
doped fiber amplifier; ESA, electrical spectrum analyzer. (B) Oscilloscope
traces of $P_{b}$ (red) and $P_{\rm{comb}}^{b}$ (blue) versus $\alpha$. (C)
Comparison of unused pump power and CE between the cases with $R=0$ (upper
panel) and $R=0.9$ (lower panel). The plot includes $P_{f}/P_{\rm{pump}}$
(solid red line), $P_{b}/P_{\rm{pump}}$ (dashed red line),
$(P_{f}+P_{b})/P_{\rm{pump}}$ (red shaded area), and CE (solid blue line) as a
function of $K$ in LLE simulation. The blue circles represent the measured CE
and the right panels shows the corresponding optical spectra of soliton
microcombs.

We explore the operation of our nanophotonic circuit by comparison to the case
with $R=0$; see Fig. 3. In both cases, the tantala PhCRs we investigate have
an intrinsic quality factor $Q_{i}=\omega_{0}/\kappa_{i}=2.7\times 10^{6}$.
Other parameters of the devices include $a=109.5$ $\mu$m, RW = 4 $\mu$m,
waveguide height = 570 nm, $A_{\rm{in}}$ = 275 nm, FSR = 200 GHz,
$D_{2}=-2\pi\times 8.5$ MHz, and $\epsilon$ = 0.9 GHz. Our setup for soliton
generation and characterization is shown in Fig. 3A; here we use an erbium-
doped fiber amplifier (EDFA) to boost the power of our tunable laser. We
characterize soliton generation by measuring at both the forward and backward
ports. To verify the solitons to be stable and free of breathing oscillations,
we photodetect a portion of the circuit output and monitor with an electrical
spectrum analyzer (ESA). Figure 3B shows the primary, real-time signals of
soliton formation in the case with a reflector and $\phi$ optimally set. As we
scan $\alpha$, we monitor $P_{b}$ (red) and $P_{\rm{comb}}^{b}$ (blue). The
soliton microcomb forms and stabilizes without an overly abrupt change in the
intraresonator power, mitigating thermal bistability. Stabilizing a soliton
with $R=0$ is marginally more sensitive to system fluctuations, such as those
of the on-chip pump power and pump detuning.

With the capability to form and measure solitons with nearly arbitrary
settings of $R$ and $\phi$, we present predictions of laser-power consumption
and CE, emphasizing how the reflector closes off excess loss in the forward
direction and enables near unit CE. We plot predictions for
$P_{f}/P_{\rm{pump}}$ (solid red line), $P_{b}/P_{\rm{pump}}$ (dashed red
line), and CE (solid blue line) as a function of $K$ in Fig. 3C. Here, the
behavior with $K$ is important because power flows in our circuit through
several channels: internal losses of the PhCR, coherent scattering to balance
the bidrectional PhCR, and the forward and backward output ports. Moreover, we
optimize $\alpha$ for maximum CE at each setting of $K$ and $P_{\rm{pump}}$.
The ratio of unused pump power, $(P_{f}+P_{b})/P_{\rm{pump}}$, indicates laser
consumption in our nanophotonic circuit as shown by the red shaded area in
Fig. 3C. The limit of $K\gg 1$, which we consider below, simplifies
understanding because the internal losses can be neglected.

Now, we directly compare the effect of the reflector. In the case with $R=0$
(upper panels of Fig. 3C), $P_{b}$ is loaded by conversion to the soliton and
$P_{f}$, and it contains about 25% of the pump for $K\gg 1$. At most 50% of
the pump is converted to the soliton, indicating clamped CE from the excess
output of double-sided emission. Indeed, our measurements of CE versus $K$
(blue circles in Fig. 3C) confirm the simulation. With the reflector set to
$R=0.9$, $P_{f}$ is substantially reduced for all values of $K$; see the lower
panels of Fig. 3C. Moreover, the backward pump inside the resonator is
enhanced because we set $\phi$ for optimal constructive interference in both
the measurements and the simulation. With appropriate $P_{\rm{pump}}$ and
$\alpha$, this one-sided nanophotonic circuit enables conversion of the
enhanced backward pump to the soliton microcomb, consuming in principle all of
the input power. Moreover, the ratio of unused pump power exhibits a rapid
decreases with $K$, reaching 0.2 at $K=2$. Even larger $K$ reduces the ratio
of unused pump power to below 0.1 and increases CE to 87%, which is ultimately
limited by the intrinsic loss of the resonator and non-unit $R$ of the
reflector. The panels labeled 1 and 2 at right in Fig. 3C show optical spectra
of soliton microcombs that correspond to the CE measurements. In both cases,
we observe an approximately constant power per mode: there are 20 comb modes
with spectral power variation below 10 dB. Particularly, in Panel 2, the
nanophotonic circuit can generate 14 comb lines with power above 0.5 mW at
$P_{\rm{pump}}$ = 33 mW. Moreover, the pump mode is 20 dB lower than the
adjacent comb lines, indicating high laser-power consumption.

Figure 4:  Practical characteristics of the nanophotonic circuit. (A)
Comparison between the measured (red circles) and simulated (solid red line)
threshold power versus $\phi$ in the nanophotonic circuit. The dashed red line
shows that of a PhCR with the same parameters except $R=0$. (B) Measured CE
(blue circles) and $P_{\rm{pump}}$ (red circles) of 11 nanophotonic circuit
devices on a single chip. (C) The optical spectrum of a 200 GHz soliton
microcomb (upper panel) and its frequency mismatch with ITU-T grid (lower
panel). (D) Frequency tuning of the DFB laser from 193.1 THz (upper panel) and
the associated spectrum generated by the nanophotonic circuit (lower panel) as
we increase the laser current from 117 mA.

Since our nanophotonic circuit generates a soliton microcomb with favorable
properties for applications, we study some of its practical characteristics,
such as operating power, device yield in fabrication, and compatibility with a
DFB laser diode as pump source; see Fig. 4. The reflector phase also
influences the threshold power for optical parametric oscillation, the
precursor of soliton formation, since an optimal setting of $\phi$ maximally
enhances the intraresonator pump. Figure 4A shows measurements (red circles)
and the LLE simulation (solid red line) of threshold power versus $\phi$ with
the reflector. In this experiment, we undercouple the resonator with $K=0.5$
to emphasize the reduced threshold power in our nanophotonic circuit while
other parameters of the devices are the same as that in Fig. 3. The threshold
power for a PhCR without a reflector (dashed red line in Fig. 4A) is a factor
of 2.9 higher than the optimal reflector circuit, since the pump power is
bidirectionally distributed not only inside the PhCR but also in the coupling
waveguide.

The capability to yield a large fraction of fabricated nanophotonic circuits
is an important consideration for applications. In particular, fabrication
imprecision can degrade performance due to the sensitivity of $\phi$ and $K$
to device geometry. We study this by testing several devices to assess
variation of CE and the $P_{\rm{pump}}$ setting used to operate the circuit;
see Fig. 4B. Here, we select 11 nanophotonic circuit devices on a single chip;
all of the devices that we fabricate are operational but the designed
parameter variations lead to reduced yield. Among the devices, CE varies
around 50% and $P_{\rm{pump}}$ needed for operation varies around 40 mW. These
values are consistent with our overall observed tolerance, and we optimize the
fabrication procedure and scanning electron microscopy inspection to maintain
a consistent $K$.

Beyond CE, microcomb applications also benefit from frequency alignment to the
ITU-T grid, the international channel standard for data communication systems.
Such alignment ensures compatibility with system components like wavelength-
division multiplexers and demultiplexers (?). We demonstrate robust frequency
alignment of the soliton microcomb generated in our nanophotonic circuit. In
Fig. 4C, we show the optical spectrum of a soliton microcomb with a precise
200 GHz mode spacing. By varying $a$ and RW of the PhCR, we finely tune both
FSR towards a target of 200 GHz and $\omega_{0}$ to align with any pre-defined
grid. We present measurements of the frequency mismatch in Fig. 4C in which
the 27 microcomb modes between 191 THz and 197 THz are aligned to the ITU-T
grid better than 2.6 GHz.

To demonstrate the compatibility of our nanophotonic circuit with chip-scale
lasers, we use a standard, butterfly-packaged DFB laser as the pump; see Fig.
4D. Such DFB lasers offer frequency adjustment through the control of laser
temperature and current; in particular, the frequency-control parameters of
our laser are -11.3 GHz/K and -0.9 GHz/mA, respectively. Hence, we tune the
DFB laser temperature to coarsely align the laser with the split mode, and
then we use the DFB laser current to sweep $\alpha$ for soliton generation.
Figure 4D illustrates the DFB laser frequency change ($\Delta f_{\rm{DFB}}$)
with current and the associated spectrum generated by the nanophotonic
circuit. We observe optical parametric oscillation followed by the evolution
to a stable soliton. This straightforward and reliable technique combined with
low required pump power and high CE of our nanophotonic circuits points the
way to a scalable, frequency-comb laser source for high-capacity
telecommunication applications.

In conclusion, we have presented a nanophotonic resonator circuit that
consists of a bidirectional Kerr resonator and a reflector on the coupling
waveguide. Coherent backscattering with nanophotonics enables arbitrary phase-
matching and high laser-power consumption in the normal dispersion PhCRs,
allowing for soliton formation at approaching unit CE. Using this nanophotonic
circuit, we generate 200 GHz flattop solitons with up to 65% CE and 97% pump
consumption. We investigate the nonlinear dynamics and steady-state fields in
the nanophotonic circuit with LLE and observe agreement between simulation and
measurement results. Moreover, we study some practical characteristics of the
nanophotonic circuit by demonstrating high yield fabrication of the high-
efficiency devices, generation of comb lines that are precisely aligned with
ITU-T frequency grids and the compatibility with a chip-scale DFB pump laser
for future integration (?, ?).

We thank Atasi Dan, Yang Li, and Daniel Slichter for reading the manuscript.
This work is a contribution of NIST and not subject to US copyright. This
research has been funded by DARPA PIPES, AFOSR FA9550-20-1-0004 Project Number
19RT1019, NSF Quantum Leap Challenge Institute Award OMA - 2016244, and NIST
on a Chip. Mention of specific companies or trade names is for scientific
communication only and does not constitute an endorsement by NIST.

## References and Notes

  * 1. A. E. Willner, S. Khaleghi, M. R. Chitgarha, O. F. Yilmaz, Journal of Lightwave Technology 32, 660 (2014).
  * 2. D. R. Carlson, et al., Science 361, 1358 (2018).
  * 3. J. Hansryd, P. Andrekson, M. Westlund, J. Li, P.-O. Hedekvist, IEEE Journal of Selected Topics in Quantum Electronics 8, 506 (2002).
  * 4. M. A. Foster, et al., Nature 441, 960 (2006).
  * 5. J. Riemensberger, et al., Nature 612, 56 (2022).
  * 6. A. Dutt, et al., Phys. Rev. Appl. 3, 044005 (2015).
  * 7. C. Reimer, et al., Science 351, 1176 (2016).
  * 8. Z. Yang, et al., Nature Communications 12, 4781 (2021).
  * 9. M. A. Guidry, D. M. Lukin, K. Y. Yang, R. Trivedi, J. Vučković, Nature Photonics 16, 52 (2022).
  * 10. M. Jahanbozorgi, et al., Optica 10, 1100 (2023).
  * 11. T. J. Kippenberg, R. Holzwarth, S. A. Diddams, Science 332, 555 (2011).
  * 12. D. J. Moss, R. Morandotti, A. L. Gaeta, M. Lipson, Nature Photonics 7, 597 (2013).
  * 13. A. Pasquazi, et al., Physics Reports 729, 1 (2018).
  * 14. T. J. Kippenberg, A. L. Gaeta, M. Lipson, M. L. Gorodetsky, Science 361, eaan8083 (2018).
  * 15. A. L. Gaeta, M. Lipson, T. J. Kippenberg, Nature Photonics 13, 158 (2019).
  * 16. X. Ji, et al., Optica 4, 619 (2017).
  * 17. X. Lu, G. Moille, A. Rao, D. A. Westly, K. Srinivasan, Optica 7, 1417 (2020).
  * 18. J. A. Black, et al., Optica 9, 1183 (2022).
  * 19. E. F. Perez, et al., Nature Communications 14, 242 (2023).
  * 20. T. Herr, et al., Nature Photonics 8, 145 (2014).
  * 21. X. Yi, Q.-F. Yang, K. Y. Yang, M.-G. Suh, K. Vahala, Optica 2, 1078 (2015).
  * 22. G. Xu, et al., Nature Communications 12, 4023 (2021).
  * 23. S.-P. Yu, et al., Nature Photonics 15, 461 (2021).
  * 24. S.-P. Yu, E. Lucas, J. Zang, S. B. Papp, Nature Communications 13, 3134 (2022).
  * 25. E. Lucas, S.-P. Yu, T. C. Briles, D. R. Carlson, S. B. Papp, Nature Photonics 17, 943 (2023).
  * 26. B. Shen, et al., Nature 582, 365 (2020).
  * 27. A. S. Voloshin, et al., Nature Communications 12, 235 (2021).
  * 28. C. Xiang, et al., Science 373, 99 (2021).
  * 29. T. C. Briles, et al., APL Photonics 6, 026102 (2021).
  * 30. A. E. Ulanov, et al., Nature Photonics (2024).
  * 31. K. Y. Yang, et al., Nature Communications 13, 7862 (2022).
  * 32. A. Rizzo, et al., Nature Photonics 17, 781 (2023).
  * 33. M. Tan, et al., Optics Communications 465, 125563 (2020).
  * 34. X. Xue, X. Zheng, B. Zhou, Nature Photonics 13, 616 (2019).
  * 35. H. Bao, et al., Nature Photonics 13, 384 (2019).
  * 36. D. C. Cole, E. S. Lamb, P. Del’Haye, S. A. Diddams, S. B. Papp, Nature Photonics 11, 671 (2017).
  * 37. Óskar B. Helgason, et al., Nature Photonics 17, 992 (2023).
  * 38. X. Xue, et al., Nature Photonics 9, 594 (2015).
  * 39. X. Xue, P.-H. Wang, Y. Xuan, M. Qi, A. M. Weiner, Laser & Photonics Reviews 11, 1600276 (2017).
  * 40. B. Y. Kim, et al., Opt. Lett. 44, 4475 (2019).
  * 41. Óskar B. Helgason, et al., Nature Photonics 15, 305 (2021).
  * 42. K. F. Lamee, D. R. Carlson, Z. L. Newman, S.-P. Yu, S. B. Papp, Optics Letters 45, 4192 (2020).
  * 43. H. Jung, et al., Optica 8, 811 (2021).
  * 44. J. A. Black, et al., Optics Letters 46, 817 (2021).
  * 45. L. A. Lugiato, R. Lefever, Phys. Rev. Lett. 58, 2209 (1987).
  * 46. H. Liu, et al., Phys. Rev. Lett. 132, 023801 (2024).
  * 47. C. Shirpurkar, et al., Optics Letters 47, 3063 (2022).


Scale-free tree network with an ultra-large diameter
Fei Maa,111 The author’s E-mail<EMAIL_ADDRESS>, Ping Wangb,c,d,222 The
corresponding author’s E-mail<EMAIL_ADDRESS>
a School of Electronics Engineering and Computer Science, Peking University,
Beijing 100871, China
b National Engineering Research Center for Software Engineering, Peking
University, Beijing, China
c School of Software and Microelectronics, Peking University, Beijing 102600,
China
d Key Laboratory of High Confidence Software Technologies (PKU), Ministry of
Education, Beijing, China

> Abstract: Scale-free networks are prevalently observed in a great variety of
> complex systems, which triggers various researches relevant to networked
> models of such type. In this work, we propose a family of growth tree
> networks $\mathcal{T}_{t}$, which turn out to be scale-free, in an iterative
> manner. As opposed to most of published tree models with scale-free feature,
> our tree networks have the power-law exponent $\gamma=1+\ln 5/\ln 2$ that is
> obviously larger than $3$. At the same time, “small-world” property can not
> be found particularly because models $\mathcal{T}_{t}$ have an ultra-large
> diameter $D_{t}$, i.e., $D_{t}\sim|\mathcal{T}_{t}|^{\ln 3/\ln 5}$ where
> $|\mathcal{T}_{t}|$ represents vertex number. In addition, we also study
> random walks on tree networks $\mathcal{T}_{t}$ and derive exact solution to
> mean hitting time $\langle\mathcal{H}_{t}\rangle$. The results suggest that
> analytic formula for quantity $\langle\mathcal{H}_{t}\rangle$ as a function
> of vertex number $\mathcal{T}_{t}$ has a power-law form, i.e.,
> $\langle\mathcal{H}_{t}\rangle\sim|\mathcal{T}_{t}|^{1+\ln 3/\ln 5}$ , a
> characteristic that is rarely encountered in the realm of scale-free tree
> networks. Lastly, we execute extensive experimental simulations, and find
> that empirical analysis is in strong agreement with theoretical results.
>
>
> Keywords: Tree, Scale-free network, Degree distribution, Diameter, Random
> walks, Mean hitting time.
>

##  1 Introduction

Nowadays, complex network, as a powerful tool, has proven useful in a great
number of applications in different fields ranging from statistic physics,
applied mathematic, computer science to chemistry, even to social science
[2]-[4]. In particular, the past two decades have seen an upsurge in the area
of complex networks. In 1998, Watts and Strogatz established the well-known
WS-model for purpose of mimicking small-world property popularly observed on
various complex systems [5]. Generally speaking, a networked model is
considered small-world if it has both smaller diameter and higher clustering
coefficient. Almost at the same time, Barabási and Albert in 1999 published
the famous BA-model in order to capture scale-free feature prevalently
encountered in a wide range of complex networks [6]. A network is called
scale-free model when it obeys power-law degree distribution. The core
ingredients in BA-model are _growth_ and _preferential attachment mechanism_.
Indeed, since then, a lot of networked models with small-world and scale-free
characteristics [7]-[11], both real-world and synthetic, have been proposed
based mainly on thought reflected by the two classic models mentioned above.
Among which, Cohen _et al_ [12] proved that scale-free networked models are
small and, sometimes, ural-small with respect to relationship between diameter
and vertex number. Afterwards, however, some theoretical scale-free models are
proved to possess a larger diameter [13]-[15]. As will be shown shortly, our
models in this work turn out to follow power-law degree distribution and,
however, have a larger diameter. To put this another way, our scale-free
models are not “small” in terms of demonstration in [12].

As known, there are various kinds of networked models proposed in the last
years, such as Apollonian Networks [16], Pseudofractal scale-free web [17],
Hierarchical networks [18] and flower graphs [19]. Tree network, as the
simplest connected networked model, has widely studied in various disciplines
[20]-[31]. For instance, the scale-free BA-tree is immediately obtained by
setting the number $m$ of edges originating from each newly added vertex in
BA-model to $1$ [6]. Clearly, the scale-free BA-tree is a stochastic model.
Similarly, Souza _et al_ [22] also proposed a family of random recursive trees
and found that greedy choice broadens the degree distribution according to the
power of choice. In [23], Gabel _et al_ discussed growth tree networks and
studied the associated structural parameters. In addition, some trees with
deterministic structure have also been studied in detail, such as,
deterministic uniform recursive tree [24], growing treelike networks [25],
Fibonacci treelike models [26], Vicsek fractal [27]-[28] and T-graph
[29]-[31]. In this paper, we will introduce a class of tree networks sharing
interesting features as will be shown later.

Fig.1. (Color online) The diagram of tree networks $\mathcal{T}_{t}$ and
$\mathcal{T}^{\star}_{t}$. The first three generations of tree
$\mathcal{T}_{t}$ are plotted in panel (a). Among them, we make use of
distinct colors to highlight vertices added at different time steps.
Similarly, we also show the first three generations of tree
$\mathcal{T}^{\star}_{t}$ in panel (b).

Often, in order to better understand the underlying structures on tree
networks, some simple yet intriguing dynamics including random walks are
considered [32]-[36]. In [27], Vicsek studied random walks on Vicsek fractal
and derived mean hitting time for random walks. Analogously, Zhang _et al_
[25] discussed mean first-passage time for random walks on growing treelike
networks using spectral technique. More generally, it is challenging for a
growth tree network with thousands of vertices to determine exact solution to
some structural parameters including mean hitting time. On the other hand, it
is of significant interest to derive closed-form solutions to topological
parameters. Roughly speaking, there are two main reasons for this. The one is
to in-depth unveil some unknown structural characteristics behind tree
networks under consideration. The other is to develop some effective methods
suitable for addressing some specific dynamical problems on networks. In this
work, we study random walks on networks proposed and precisely obtain what we
are interested in by proposing a series of simple combinatorial methods
instead of the commonly used spectral technique.

The rest of this paper is organized as follows. Section II aims at introducing
the concrete procedure for generating our tree network $\mathcal{T}_{t}$. At
the same time, vertex number is obtained in an iterative manner. Next, we
estimate two significant structural parameters, i.e., degree distribution and
diameter, in Section III. The former suggests that tree network
$\mathcal{T}_{t}$ follows power-law degree with exponent larger than $3$. The
latter indicates that “small-world” property is not observed on network
$\mathcal{T}_{t}$. In Section IV, we study random walks on tree network
$\mathcal{T}_{t}$, and derive closed-form solution of mean hitting time. The
corresponding analytic formula for mean hitting time as a function of vertex
number exhibits a power-law form. Note that throughout the preceding two
sections, we also in-depth discuss the adjoint tree model of tree network
$\mathcal{T}_{t}$ which is viewed as $\mathcal{T}^{\star}_{t}$. Finally, we
close this work by drawing the conclusion in Section V.

##  2 Construction

In this section, we propose an algorithm $\Phi$ for creating tree network
$\mathcal{T}_{t}$, and then derive some fundamental parameters. Part of them
are helpful to analyze other structural quantities as will be shown later. At
the meantime, the corresponding adjoint tree network $\mathcal{T}^{\star}_{t}$
is also established. Note that the exact solutions of some structural
parameters relevant to tree $\mathcal{T}^{\star}_{t}$ are only given without
detailed derivations. This is because they can be obtained in a similar manner
to that for tree network $\mathcal{T}_{t}$.

###  2.1 Algorithm $\Phi$

First of all, it should be mentioned that an arbitrary tree $\mathcal{T}$ can
be selected as input to the algorithm $\Phi$ proposed below. Nonetheless, as
tried in the rich literature [25, 27, 30], we only consider a specific case,
i.e., a single edge serving as the seed of algorithm $\Phi$, when studying
degree distribution. In other cases, we still consider an arbitrary tree as
seed.

At $t=0$, the seminal model is an arbitrary tree $\mathcal{T}$, referred to as
$\mathcal{T}_{0}$. For our purpose, we denote by $|\mathcal{T}_{0}|$ vertex
number.

At $t=1$, the young tree $\mathcal{T}_{1}$ is obtained from parent tree
$\mathcal{T}_{0}$ by performing the following two operations:

  * •

We first insert one new vertex into each edge incident with each vertex $v$
belonging to tree $\mathcal{T}_{0}$.

  * •

And then, we connect $k_{v}(0)$ new vertices as leaves to each vertex $v$ of
tree $\mathcal{T}_{0}$. Here, symbol $k_{v}(0)$ represents the degree of
vertex $v$.

_Remark_ —Such a manipulation is easy to extend to contain many more
complicated operations. Yet, our goal is to introduce this idea and thus some
more intriguing extensions are left as investigations in future research.

At $t\geq 2$, the next tree $\mathcal{T}_{t}$ is generated based on the
previous tree $\mathcal{T}_{t-1}$ using the same operation mentioned above.

To make our description more concrete, an illustrative example is shown in
panel (a) of Fig.1 where a single edge is chose as a seed of algorithm $\Phi$.

###  2.2 Vertex number

Clearly, our tree network $\mathcal{T}_{t}$ is also built upon an iterative
method, a generative approach that has been frequently used to create various
networked models in the past years [9, 13, 17, 18]. This enables calculations
of some fundamental topological parameters on network created, such as, vertex
number. Below, we provide a simple derivation of vertex number
$|\mathcal{T}_{t}|$.

From algorithm $\Phi$, we can see that in fact, there are $2k_{v}(t-1)$ new
vertices created for each vertex $v$ of tree $\mathcal{T}_{t-1}$ at time step
$t$. This leads to the following relation

$|\mathcal{T}_{t}|=\sum_{v\in\mathcal{T}_{t-1}}(1+2k_{v}(t-1)).$

With a trivial fact that the summation of degrees of all vertices in a tree is
equal to two times (vertex number minus one), the exact solution to vertex
number $|\mathcal{T}_{t}|$ is derived by calculating a recurrence and then is
given by

$|\mathcal{T}_{t}|=5^{t}(|\mathcal{T}_{0}|-1)+1.$ (1)

Besides that, we want to make more detailed demonstration about algorithm
$\Phi$ in order to enable the next calculations. Specifically, in practice,
there are two classes of new vertices added at each time step. As a case
study, it is easy to see that $k_{v}(t-1)$ new vertices are inserted on all
the existing edges incident with vertex $v$ of $\mathcal{T}_{t-1}$ with each
being in each edge. Thus, each of vertices of this type has initial degree
$2$. On the contrary, the other $k_{v}(t-1)$ new vertices are directly
connected to vertex $v$ itself, which means that each of them has initial
degree $1$. For brevity, the former vertices associated with vertex $v$ are
grouped into set $A_{v}(t)$, and, similarly, notation $B_{v}(t)$ is used to
indicate the set of the latter vertices. To put this another way, at time step
$t$, tree network $\mathcal{T}_{t}$ has three types of vertex sets, i.e.,
$\mathcal{T}_{t-1}$, $\bigcup_{v\in\mathcal{T}_{t-1}}A_{v}(t)$ and
$\bigcup_{v\in\mathcal{T}_{t-1}}B_{v}(t)$. To make further progress, we define
$A(t)=\bigcup_{v\in\mathcal{T}_{t-1}}A_{v}(t)$ and
$B(t)=\bigcup_{v\in\mathcal{T}_{t-1}}B_{v}(t)$.

So far, we finish the construction of tree network $\mathcal{T}_{t}$. By the
way, if we adjust algorithm $\Phi$ to let $k_{v}(t-1)$ new vertices in set
$A_{v}(t)$ not be placed on edges incident with vertex $v$ but also connect to
vertex $v$ itself directly at each time step $t$, then we obtain another tree
network $\mathcal{T}^{\star}_{t}$. Panel (b) of Fig.1 shows the first three
generations of tree network $\mathcal{T}^{\star}_{t}$. It is obvious to
understand that tree network $\mathcal{T}^{\star}_{t}$ shares an identical
vertex number with tree network $\mathcal{T}_{t}$ when a common tree is
selected as their seeds, i.e.,

$|\mathcal{T}^{\star}_{t}|=5^{t}(|\mathcal{T}_{0}|-1)+1.$ (2)

Nonetheless, the both tree networks still exhibit significant difference on
underlying structure when we discuss some other structural parameters as will
be shown shortly.

##  3 Two significant topological parameters

From now on, we will focus on two significant structural parameters, degree
distribution and diameter, on tree network $\mathcal{T}_{t}$. One of important
reasons for this is that the two parameters have been widely used to study
various networked models in the literature [2]-[6]. It is using degree
distribution that Barabási _et al_ revealed the scale-free feature of a wide
range of networks [6]. With diameter, a great number of networks turn out to
be small-world [5]. More generally, determining structural parameters on
networked models is useful to understand underlying structure on networks. As
will see in the following, our tree network $\mathcal{T}_{t}$ indeed shows
some characteristics associated with both degree distribution and diameter
that can not be found in some previous tree models [24]-[31].

###  3.1 Degree distribution

Degree distribution, as an important index, has been popularly adopted in the
study of complex networks [2]. Informally, ones always use $P(k)$ to define
the probability of choosing at random uniformly a degree $k$ vertex in a given
network $\mathcal{G}(\mathcal{V},\mathcal{E})$. In the thermodynamic limit,
probability $P(k)$ is referred to as the ratio of the total number $N_{k}$ of
vertices with degree $k$ and vertex number $|\mathcal{V}|$. According to the
fact that vertex degrees of network $\mathcal{G}(\mathcal{V},\mathcal{E})$ are
discrete values in form, the accumulative degree distribution $P_{cum}(k)$ is
often exploited, which is by definition written as

$P_{cum}(k)=\frac{\sum_{k^{\prime}\geq k}N_{k^{\prime}}}{|\mathcal{V}|}.$

Clearly, it suffices to rank all the vertices of tree network
$\mathcal{T}_{t}$ with respect to vertex degree in order to derive
accumulative degree distribution $P_{cum}(k)$. Now, let us keep this in mind
and then recall the concrete procedure in algorithm $\Phi$. For one vertex $v$
added at time step $t_{i}$ into tree $\mathcal{T}_{t}$, we can see that its
degree $k_{v}(t)$ might take two distinct values. If vertex $v$ is in set
$A(t_{i})$, then the corresponding degree $k_{v}(t)$ is equal to
$2^{t-t_{i}+1}$. On the other hand, degree $k_{v}(t)$ is equivalent to
$2^{t-t_{i}}$ when vertex $v$ belongs to set $B(t_{i})$. Based on this, we can
obtain the following expression

$P_{cum}(k)=\frac{2+\sum_{j=1}^{t_{i}}|A(j)|+\sum_{j=1}^{t_{i}-1}|B(j)|}{5^{t}+1}$
(3)

where $|A|$ indicates the cardinality of a set $A$ and we have made use of a
single edge as seed. Substituting $k=2^{t-t_{i}+1}$ into Eq.(12) yields

$P_{cum}(k)\sim k^{-\gamma},\quad\gamma=\frac{\ln 5}{\ln 2}.$

Performing derivative on both hand sides in the above expression produces

$P(k)\sim k^{-(1+\gamma)},$ (4)

which suggests that our tree $\mathcal{T}_{t}$ follows power-law degree
distribution. In other words, scale-free feature is observed on tree
$\mathcal{T}_{t}$ in the limit of large graph size. It is worth noting that
power-law exponent $1+\gamma$ is certainly larger than $3$. This is a
phenomenon that is seldom reported in published tree networks with scale-free
feature. In essence, it is also confirmed in terms of “preferential
attachment” in the following observation.

Analogously, we can without difficulty to derive degree distribution of tree
$\mathcal{T}^{\star}_{t}$, which is showed

$P^{\star}(k)\sim k^{-(1+\gamma^{\star})},\quad\gamma^{\star}=\frac{\ln 5}{\ln
3}.$ (5)

This means that tree $\mathcal{T}^{\star}_{t}$ obeys power-law degree
distribution as well. However, its power-law exponent is smaller than $3$.

Apparently, degree distribution implies that tree $\mathcal{T}^{\star}_{t}$
has a different underlying structure from tree $\mathcal{T}_{t}$. Yet, it
should be mentioned that in the large graph size limit, they have an identical
average degree regardless of power-law exponent. In addition, the both types
of generative ways to create scale-free tree networks $\mathcal{T}_{t}$ and
$\mathcal{T}^{\star}_{t}$ are within the category of “preferential attachment
mechanism” introduced by Barabási _et al_ [6]. In the past two decades, most
papers focus on the latter generative approach to generating tree
$\mathcal{T}^{\star}_{t}$ [26, 37]. As a result, almost all tree networks with
scale-free feature are always proved to have an exponent no more than $3$.
Therefore, this work provides an effective way to build up scale-free tree
networks with intriguing features, such as power-law exponent larger than $3$.

###  3.2 Diameter

In the jargon of graph theory, diameter of network
$\mathcal{G}(\mathcal{V},\mathcal{E})$, denoted by $D$, is the distance of
longest shortest path between all possible pairs of vertices. As one of most
fundamental topological parameters, diameter has been proved useful in a
variety of applications [38]-[40]. For example, it is in general utilized to
measure information delay in communication network. As known, small-world
phenomena on complex networks are certified by means of diameter [5]. Next,
let us pay attention to discussion about diameter of tree network
$\mathcal{T}_{t}$.

As shown in algorithm $\Phi$, two new vertices are inserted into each pre-
existing edge at each time step. This results in the following relation
satisfied by diameter $D_{t}$ of tree $\mathcal{T}_{t}$, namely,

$D_{t}=2+3D_{t-1}.$

With initial condition $D_{0}$ that is diameter of the seed $\mathcal{T}_{0}$,
it is straightforward to write

$D_{t}=3^{t}(D_{0}+1)-1.$ (6)

Obviously, unlike statement in [12], diameter $D_{t}$ of tree network
$\mathcal{T}_{t}$ is much more than $\ln|\mathcal{T}_{t}|$ and follows
$D_{t}\sim|\mathcal{T}_{t}|^{\ln 3/\ln 5}$. This means that our tree network
$\mathcal{T}_{t}$ has no small-world property.

As above, we also obtain exact solution of diameter $D^{\star}_{t}$ of tree
network $\mathcal{T}^{\star}_{t}$ in an iterative fashion and write

$D^{\star}_{t}=D_{0}+2t$ (7)

in which $D_{0}$ has the same meaning as mentioned above. It is clear to see
that tree $\mathcal{T}^{\star}_{t}$ meets the criteria in [12], i.e.,
$D^{\star}_{t}\sim\ln|\mathcal{T}^{\star}_{t}|$. In some published papers
[25], trees of such type are considered as models with small-world property.
Strictly speaking, an arbitrary tree network can not be small-world mainly due
to the corresponding concept, i.e., small diameter and higher clustering
coefficient. Tree network has zero clustering coefficient.

From Eqs.(6) and (7), we can also state that the both types of generative ways
have significantly different impacts on structure of tree networks.
Intuitively, the way to create tree $\mathcal{T}_{t}$ enlarges the underlying
structure of network more remarkably than that for creating tree
$\mathcal{T}^{\star}_{t}$. Therefore, it is a useful manner in which one may
build up some networked models having a quite large diameter. Note also that a
similar operation, which makes the growth of diameter of network under
consideration more rapid, has been used to construct many other networked
models of great interest [13]-[15], including scale-free networks with both
density feature and large diameter [41].

##  4 Random walks on tree network

Random walks, as the simplest representative of discrete-time unbiased
diffusion phenomena, have been receiving increasing attention from various
fields ranging from applied mathematics, statistic physics, computer science
to bio-chemistry [29, 33, 34, 35, 42, 43, 44], and so forth. A walker
performing random walks on network $\mathcal{G}(\mathcal{V},\mathcal{E})$
moves with a probability $P_{u\rightarrow v}=1/k_{u}$ from his current
position $u$ to each candidate $v$ in the neighboring set $\Omega_{u}$ in a
step. Such a dynamic process is in general encoded by the corresponding
transition matrix $\mathbf{P}=\mathbf{D}^{-1}\mathbf{A}$ of network
$\mathcal{G}(\mathcal{V},\mathcal{E})$ where $\mathbf{A}$ is adjacency matrix
whose entry $a_{uv}$ is given by

$a_{uv}=\left\\{\begin{aligned} &1,\quad\text{vertex $u$ is connected to $v$
by an edge}\\\ &0,\quad\text{otherwise}.\end{aligned}\right.$

and $\mathbf{D}$ is the diagonal matrix that is defined as below: the $u$-th
diagonal entry is $k_{u}$, while all non-diagonal entries are zero, i.e.,
$\mathbf{D}=\text{diag}[k_{1},k_{2},\dots,k_{|\mathcal{V}|}]$. Based on this,
the significant structural parameter, called hitting time that is denoted by
$\mathcal{H}_{u\rightarrow v}$, is written as

$\mathcal{H}_{u\rightarrow v}=\left.\dfrac{d}{dz}\mathcal{P}_{u\rightarrow
v}(z)\right|_{z=1},$ (8)

where $\mathcal{P}_{u\rightarrow v}(z)$ is the probability generating function
that is referred to as

$\mathcal{P}_{u\rightarrow v}(z)=\sum_{t=0}^{\infty}P_{u\rightarrow
v}(t)z^{t}.$

Here, we have made use of the probability $P_{u\rightarrow v}(t)$ that a
walker starting out from source vertex $u$ first reaches to his destination
vertex $v$ after $t$ steps. For network $\mathcal{G}(\mathcal{V},\mathcal{E})$
as a whole, mean hitting time, denoted by $\langle\mathcal{H}\rangle$, is
viewed as

$\langle\mathcal{H}\rangle=\frac{1}{|\mathcal{V}|(|\mathcal{V}|-1)}\sum_{u,v}\mathcal{H}_{u\rightarrow
v}.$ (9)

Note also that the quantity $\mathcal{H}_{u\rightarrow v}$ is not necessarily
equal to $\mathcal{H}_{v\rightarrow u}$ in general. As known, however,
equality between these two quantities holds for tree network. What is even
more interesting is that for a tree $\mathcal{T}$, the summation over the both
quantities is closely connected to their own distance $d_{uv}$ in the
following term

$\mathcal{H}_{u\rightarrow v}+\mathcal{H}_{v\rightarrow
u}=2d_{uv}(|\mathcal{T}|-1).$ (10)

What’s more, a more general version associated with Eq.(10) has been reported
in [45, 46]. In fact, such a summation in left-hand side of Eq.(10) is the so-
called commute time between vertices $u$ and $v$. And then, mean hitting time
$\langle\mathcal{H}_{\mathcal{T}}\rangle$ on tree network [35] can be
reorganized as

$2\mathcal{W}_{\mathcal{T}}=|\mathcal{T}|\langle\mathcal{H}_{\mathcal{T}}\rangle,$
(11)

where $\mathcal{W}_{\mathcal{T}}$ is Wiener index of tree $\mathcal{T}$, and
is by definition expressed as

$\mathcal{W}_{\mathcal{T}}=\frac{1}{2}\sum_{u,v\in\mathcal{T}}d_{uv}.$

Clearly, determining the exact solution to mean hitting time on tree network
is equivalently transformed into the problem of calculating the corresponding
Wiener index. This enables calculation of mean hitting time on tree when it is
more convenient to deriving the latter, and vice verse. As will see, the
following derivation of mean hitting time on our tree network
$\mathcal{T}_{t}$ is carried out based on Eq.(11).

###  4.1 Mean hitting time

In this section, we will precisely determine the closed-form solution to mean
hitting time on tree networks $\mathcal{T}_{t}$ and $\mathcal{T}^{\star}_{t}$.
Specifically, a rigorous derivation for mean hitting time
$\langle\mathcal{H}_{t}\rangle$ for random walks on tree $\mathcal{T}_{t}$ is
demonstrated in stages. However, we immediately give an approximate solution
to $\langle\mathcal{H}^{\star}_{t}\rangle$ on tree $\mathcal{T}^{\star}_{t}$
due to a similar discussion to that for quantity
$\langle\mathcal{H}_{t}\rangle$.

Beginning with our calculations, let us recall the concrete procedure for
developing scale-free tree network $\mathcal{T}_{t}$. At time step $t$, tree
$\mathcal{T}_{t}$ consists of three classes of vertex sets, i.e.,
$\mathcal{T}_{t-1}$, $A(t)=\bigcup_{v\in\mathcal{T}_{t-1}}A_{v}(t)$ and
$B(t)=\bigcup_{v\in\mathcal{T}_{t-1}}B_{v}(t)$. It is natural to see that by
definition, Wiener index $\mathcal{W}_{t}$ of tree $\mathcal{T}_{t}$ is
obtained if we can determine distances of all possible pairs of vertices from
the three sets. Roughly speaking, there are nine distinct combinatorial
manners between vertices according to different classifications of vertices of
tree $\mathcal{T}_{t}$. For brevity and our purpose, given a vertex $v$ of
tree $\mathcal{T}_{t-1}$, the corresponding each vertex in set $A_{v}(t)$ is
labeled $v^{1}_{i}$ and, similarly, each vertex in set $B_{v}(t)$ is labeled
$v^{2}_{i}$. Now, we are ready to estimate quantity
$\langle\mathcal{H}_{t}\rangle$ on tree $\mathcal{T}_{t}$.

_Stage 1._ Due to the fact that three are only two new vertices inserted on
each edge $uv$ of tree $\mathcal{T}_{t-1}$, we can without difficulty to
obtain

$\mathcal{W}_{t}(1)=\frac{1}{2}\sum_{u,v\in\mathcal{T}_{t-1}}d^{\dagger}_{uv}=3\mathcal{W}_{t-1}.$
(12)

Here, $d^{\dagger}_{uv}$ represents distance between vertices $u$ and $v$ when
considering tree $\mathcal{T}_{t}$. It is worth noticing that, hereafter, we
make use of the superscribe $\dagger$ to indicate distance between vertex pair
on tree $\mathcal{T}_{t}$.

_Stage 2._ For each vertex $v$ in tree $\mathcal{T}_{t-1}$, there are in
practice $2k_{v}(t-1)$ vertices added into tree $\mathcal{T}_{t}$ where
$k_{v}(t-1)$ is the degree of vertex $v$ in tree $\mathcal{T}_{t-1}$. Half of
them are grouped into set $A_{v}(t)$. The other is in set $B_{v}(t)$. Using
notations defined above, it is easy to see the following two equations

$\displaystyle W_{t}(v;1)=\sum_{v^{1}_{i}\in
A_{v}(t)}d^{\dagger}_{vv^{1}_{i}}=k_{v}(t-1),$ (13a)
$W_{t}(v;2)=\sum_{v^{2}_{i}\in B_{v}(t)}d^{\dagger}_{vv^{2}_{i}}=k_{v}(t-1).$
(13b)

In a word, we have

$\begin{aligned}
\mathcal{W}_{t}(2)&=\sum_{v\in\mathcal{T}_{t-1}}\sum_{v^{1}_{i}\in
A_{v}(t)}d^{\dagger}_{vv^{1}_{i}}+\sum_{v\in\mathcal{T}_{t-1}}\sum_{v^{2}_{i}\in
A_{v}(t)}d^{\dagger}_{vv^{2}_{i}}\\\
&=\sum_{v\in\mathcal{T}_{t-1}}\left(W_{t}(v;1)+W_{t}(v;2)\right)\\\
&=2\sum_{v\in\mathcal{T}_{t-1}}k_{v}(t-1)\\\
&=4(|\mathcal{T}_{t-1}|-1)\end{aligned}.$ (14)

_Stage 3._ Following the above-mentioned stage, let us still concentrate on
local structure associated with each vertex $v$ in tree $\mathcal{T}_{t-1}$,
i.e., $A_{v}(t)$ and $B_{v}(t)$. And then, it is easy to find

$\displaystyle\begin{aligned}
W_{t}(A_{v}(t))&=\frac{1}{2}\sum_{v^{1}_{i},v^{1}_{j}\in
A_{v}(t)}d^{\dagger}_{v^{1}_{i}v^{1}_{j}}\\\
&=k_{v}(t-1)(k_{v}(t-1)-1)\end{aligned},$ (15a) $\begin{aligned}
W_{t}(B_{v}(t))&=\frac{1}{2}\sum_{v^{2}_{i},v^{2}_{j}\in
B_{v}(t)}d^{\dagger}_{v^{2}_{i}v^{2}_{j}}\\\
&=k_{v}(t-1)(k_{v}(t-1)-1)\end{aligned}.$ (15b)

Then, we write

$\displaystyle\begin{aligned}
\mathcal{W}_{t}(3)&=\sum_{v\in\mathcal{T}_{t-1}}W_{t}(A_{v}(t))\\\
&=\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1)-2(|\mathcal{T}_{t-1}|-1)\end{aligned},$
(16a) $\begin{aligned}
\mathcal{W}_{t}(3)&=\sum_{v\in\mathcal{T}_{t-1}}W_{t}(B_{v}(t))\\\
&=\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1)-2(|\mathcal{T}_{t-1}|-1)\end{aligned}.$
(16b)

_Stage 4._ Now, we attempt to study distance $d^{\dagger}_{uv^{1}_{i}}$
between vertex $u$ from tree $\mathcal{T}_{t-1}$ and vertex $v^{1}_{i}$ in set
$A_{v}(t)$. Note that vertex $u$ is distinct with $v$. In another word, vertex
$v^{1}_{i}$ does not belong to set $A_{u}(t)$. To this end, we need to
consider a trivial fact that there is a unique path between an arbitrarily
given pair of vertices, say $u$ and $v$, in tree. For convenience, we denote
by $\mathcal{P}_{uv}$ the path joining vertex $u$ to $v$. According to
underlying structure of tree $\mathcal{T}_{t}$, path
$\mathcal{P}_{uv^{1}_{i}}$ is constituted in two manners. The one is that path
$\mathcal{P}_{uv^{1}_{i}}$ consists of two segments $\mathcal{P}_{uv}$ and
$\mathcal{P}_{vv^{1}_{i}}$. The number of paths of such type is
$(k_{v}(t-1)-1)$ in total. The other is that path $\mathcal{P}_{uv^{1}_{i}}$
is obtained from path $\mathcal{P}_{uv}$ by deleting segment
$\mathcal{P}_{v^{1}_{i}v}$. This case is encountered in tree $\mathcal{T}_{t}$
only once. Intuitively, given a pair of vertices $u$ and $v$ in tree
$\mathcal{T}_{t-1}$, it is straightforward to find the following expression

$\begin{aligned} W_{t}(u;A_{v}(t))&=\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}\sum_{v^{1}_{i}\in A_{v}(t)}d^{\dagger}_{uv^{1}_{i}}\\\
&=\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}\left[k_{v}(t-1)d^{\dagger}_{uv}+(k_{v}(t-1)-2)\right]\end{aligned}.$
(17)

Summing on both hand sides of Eq.(17) with respect to all vertices $u$ in tree
$\mathcal{T}_{t-1}$ yields

$\begin{aligned} W_{t}(4)&=\sum_{u\in\mathcal{T}_{t-1}}W_{t}(u;A_{v}(t))\\\
&=\sum_{u\in\mathcal{T}_{t-1}}\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}k_{v}(t-1)d^{\dagger}_{uv}\\\
&\quad+\sum_{u\in\mathcal{T}_{t-1}}\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}(k_{v}(t-1)-2)\end{aligned}.$ (18)

However, it is not easy to derive an exact solution for quantity $W_{t}(4)$ in
terms of Eq.(18). To get around this issue, we introduce another calculational
technique as below. Specifically, we take an arbitrary path $\mathcal{P}_{uv}$
in tree $\mathcal{T}_{t}$ as an illustrative example. If the length of path
$\mathcal{P}_{uv}$ is equal to $3$, we can express
$\mathcal{P}_{uv}:=uu^{1}_{i}v^{1}_{j}v$. In this case, it is easy to obtain

$d^{\dagger}_{uv^{1}_{j}}+d^{\dagger}_{vu^{1}_{i}}=4.$ (19)

Due to nature of tree, there are number $(|\mathcal{T}_{t-1}|-1)$ of paths of
this kind in total. Besides that, if path $\mathcal{P}_{uv}$ has length larger
than $3$, then, without loss of generality, path $\mathcal{P}_{uv}$ is
expressed as

$\mathcal{P}_{uv}:=uu^{1}_{i}s^{1}_{j}s\dots ww^{1}_{l}v^{1}_{n}v.$

This leads to an equality

$d^{\dagger}_{uv^{1}_{n}}+d^{\dagger}_{uw^{1}_{l}}+d^{\dagger}_{vu^{1}_{i}}+d^{\dagger}_{vs^{1}_{j}}=4d^{\dagger}_{uv}-6.$
(20)

Armed with the results in Eq.(19) and (20), quantity $W_{t}(4)$ is rewritten
as

$\begin{aligned}
W_{t}(4)&=4[\mathcal{W}_{t}(1)-3(|\mathcal{T}_{t-1}|-1)]-6\left(\begin{array}[]{c}|\mathcal{T}_{t-1}|-1\\\
2\\\ \end{array}\right)\\\ &\quad+4(|\mathcal{T}_{t-1}|-1)\\\
&=4\mathcal{W}_{t}(1)-(|\mathcal{T}_{t-1}|-1)(3|\mathcal{T}_{t-1}|+2)\end{aligned}.$
(21)

Taking into account a fact that the both terms in the third line are easy to
answer, the closed-form solution of $W_{t}(4)$ is more conveniently calculated
in view of Eq.(21) than Eq.(18). It is clear to the eye that such a
transformation from Eq.(18) to Eq.(21) is quite effective. At the same time,
it should be mentioned that the light shed by Eq.(21) will be proved helpful
in the following derivations.

_Stage 5._ Let us divert our insight into discussion about distance
$d^{\dagger}_{u^{1}_{i}v^{1}_{j}}$ between two vertices, say $u^{1}_{i}$ and
$v^{1}_{j}$, from sets $A_{u}(t)$ and $A_{v}(t)$. It is worth noticing that
subscript $u$ can be the same as $v$. In other words, the result in Eq.(16a)
is covered by the following formula. By definition, given a couple of distinct
sets $A_{u}(t)$ and $A_{v}(t)$, the summation
$\mathcal{W}_{t}(A_{u}(t);A_{v}(t))$ of distances
$d^{\dagger}_{u^{1}_{i}v^{1}_{j}}$ of all possible vertex pairs is given by

$\displaystyle\mathcal{W}_{t}(A_{u}(t);A_{v}(t))$
$\displaystyle=k_{u}(t-1)k_{v}(t-1)d^{\dagger}_{uv}-2$ (22)
$\displaystyle\quad+2(k_{u}(t-1)-1)(k_{v}(t-1)-1)$

where $A_{u}(t)$ is distinct with $A_{v}(t)$.

To make further progress, from Eqs.(16a) and (22), we can write

$\begin{aligned}
W_{t}(5)&=\frac{1}{2}\sum_{u,v\in\mathcal{T}_{t-1}}\mathcal{W}_{t}(A_{u}(t);A_{v}(t))+\sum_{v\in\mathcal{T}_{t-1}}W_{t}(A_{v}(t))\\\
&=\frac{1}{2}\sum_{u,v\in\mathcal{T}_{t-1}}k_{u}(t-1)k_{v}(t-1)d^{\dagger}_{uv}-\left(\begin{array}[]{c}|\mathcal{T}_{t-1}|\\\
2\\\ \end{array}\right)\\\
&\quad+\sum_{u,v\in\mathcal{T}_{t-1}}(k_{u}(t-1)-1)(k_{v}(t-1)-1)\\\
&\quad+\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1)-2(|\mathcal{T}_{t-1}|-1)\end{aligned}.$
(23)

In general, determining solution to quantity $W_{t}(5)$ according to Eq.(23)
is a tough task. To address this issue, we still attempt to take advantage of
transformational method. To be more concrete, the detailed calculations are
deduced based on a similar manner to that used in stage 4. That is to say, we
continue to focus on path $\mathcal{P}_{uv}$ and also make use of some
notations mentioned above. There are two sub-cases contained in the following
discussions. For path $\mathcal{P}_{uv}$ with length $3$, distance
$d^{\dagger}_{u^{1}_{i}v^{1}_{j}}$ is calculated equal to $1$. Additionally,
if the length of path $\mathcal{P}_{uv}$ is greater than $3$, we have an
equality like that in Eq.(20) that is given by

$d^{\dagger}_{u^{1}_{i}v^{1}_{n}}+d^{\dagger}_{u^{1}_{i}w^{1}_{l}}+d^{\dagger}_{s^{1}_{j}v^{1}_{n}}+d^{\dagger}_{s^{1}_{j}w^{1}_{l}}=4d^{\dagger}_{uv}-12.$
(24)

Analogously, we can re-express quantity $W_{t}(5)$ in terms of the next
equation

$\begin{aligned}
W_{t}(5)&=4[\mathcal{W}_{t}(1)-3(|\mathcal{T}_{t-1}|-1)]-12\left(\begin{array}[]{c}|\mathcal{T}_{t-1}|-1\\\
2\\\ \end{array}\right)\\\ &\quad+(|\mathcal{T}_{t-1}|-1)\\\
&=4\mathcal{W}_{t}(1)-(|\mathcal{T}_{t-1}|-1)(6|\mathcal{T}_{t-1}|-1)\end{aligned}.$
(25)

_Stage 6._ Here, we pay more attention on determination of distance
$d^{\dagger}_{uv^{2}_{i}}$. Note that vertex $u$ is different from $v$. Along
the same line in stage 4, the summation $W_{t}(u;B_{v}(t))$ over distances
between vertex $u$ and each vertex $v^{2}_{i}$ in set $B_{v}(t)$ is calculated
to yield

$\begin{aligned} W_{t}(u;B_{v}(t))&=\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}\sum_{v^{2}_{i}\in B_{v}(t)}d^{\dagger}_{uv^{2}_{i}}\\\
&=\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}\left[k_{v}(t-1)d^{\dagger}_{uv}+k_{v}(t-1)\right]\end{aligned}.$
(26)

In order to simplify the next derivation, we now build up a bridge between
quantities $W_{t}(u;B_{v}(t))$ and $W_{t}(u;A_{v}(t))$, which is

$W_{t}(u;B_{v}(t))=W_{t}(u;A_{v}(t))+\sum_{v(\neq u)\in\mathcal{T}_{t-1}}2.$
(27)

And then, we have

$\begin{aligned} W_{t}(6)&=\sum_{u\in\mathcal{T}_{t-1}}W_{t}(u;B_{v}(t))\\\
&=\sum_{u\in\mathcal{T}_{t-1}}\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}k_{v}(t-1)d^{\dagger}_{uv}\\\
&\quad+\sum_{u\in\mathcal{T}_{t-1}}\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}k_{v}(t-1)\\\
&=W_{t}(4)+\sum_{u\in\mathcal{T}_{t-1}}\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}2\\\
&=W_{t}(4)+2|\mathcal{T}_{t-1}|(|\mathcal{T}_{t-1}|-1)\end{aligned}.$ (28)

With Eq.(21), we can obtain

$W_{t}(6)=4\mathcal{W}_{t}(1)-(|\mathcal{T}_{t-1}|-1)(|\mathcal{T}_{t-1}|+2).$
(29)

_Stage 7._ For distance $d^{\dagger}_{u^{2}_{i}v^{2}_{j}}$ between two
vertices, say $u^{2}_{i}$ and $v^{2}_{j}$, from distinct sets $B_{u}(t)$ and
$B_{v}(t)$, using a similar deduction to that used to develop quantity
$\mathcal{W}_{t}(A_{u}(t);A_{v}(t))$ in stage 5 produces

$\begin{aligned}
\mathcal{W}_{t}(B_{u}(t);B_{v}(t))&=k_{u}(t-1)k_{v}(t-1)d^{\dagger}_{uv}\\\
&\quad+2k_{u}(t-1)k_{v}(t-1)\end{aligned}.$ (30)

As above, we also find a close connection of quantity
$\mathcal{W}_{t}(B_{u}(t);B_{v}(t))$ to $\mathcal{W}_{t}(A_{u}(t);A_{v}(t))$,
namely,

$\begin{aligned}
\mathcal{W}_{t}(B_{u}(t);B_{v}(t))&=\mathcal{W}_{t}(A_{u}(t);A_{v}(t))\\\
&\quad+2(k_{u}(t-1)+k_{v}(t-1))\end{aligned}.$ (31)

After that, with the help of Eqs.(25) and (31), we obtain

$\begin{aligned}
\mathcal{W}_{t}(7)&=\frac{1}{2}\sum_{u,v\in\mathcal{T}_{t-1}}\mathcal{W}_{t}(B_{u}(t);B_{v}(t))+\sum_{v\in\mathcal{T}_{t-1}}W_{t}(B_{v}(t))=\mathcal{W}_{t}(5)+\sum_{u\in\mathcal{T}_{t-1}}\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}(k_{u}(t-1)+k_{v}(t-1))\\\
&=\mathcal{W}_{t}(5)+\sum_{u\in\mathcal{T}_{t-1}}\left[(|\mathcal{T}_{t-1}|-2)k_{u}(t-1)+2(|\mathcal{T}_{t-1}|-1)\right]\\\
&=4\mathcal{W}_{t}(1)-(|\mathcal{T}_{t-1}|-1)(2|\mathcal{T}_{t-1}|+3)\end{aligned}.$
(32)

Notice also that Eq.(16b) has been used.

_Stage 8._ The last part of our tasks is to measure distance
$d^{\dagger}_{u^{1}_{i}v^{2}_{j}}$ of vertices $u^{1}_{i}$ and $v^{2}_{j}$.
This may be finished in the following two cases. When the both vertices are
connected to an identical vertex $v$, it is clear to find an equality

$\mathcal{W}_{t}(A_{v}(t);B_{v}(t))=\sum_{v^{1}_{i}\in
A_{v}(t)}\sum_{v^{2}_{i}\in
B_{v}(t)}d^{\dagger}_{v^{1}_{i}v^{2}_{j}}=2k^{2}_{v}(t-1).$ (33)

And then, we establish

$\mathcal{W}_{t}(8)=\sum_{v\in\mathcal{T}_{t-1}}\mathcal{W}_{t}(A_{v}(t);B_{v}(t))=2\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1).$
(34)

On the other hand, if vertices $u^{1}_{i}$ and $v^{2}_{j}$ are from distinct
sets $A_{u}(t)$ and $B_{v}(t)$, respectively, then we build up

$\begin{aligned} \mathcal{W}_{t}(A_{u}(t);B_{v}(t))&=\sum_{u^{1}_{i}\in
A_{u}(t)}\sum_{v^{2}_{i}\in B_{v}(t)}d^{\dagger}_{u^{1}_{i}v^{2}_{j}}\\\
&=k_{u}(t-1)k_{v}(t-1)d^{\dagger}_{uv}\\\ &\quad+2k_{u}(t-1)(k_{v}(t-1)-1)\\\
&=\mathcal{W}_{t}(B_{u}(t);B_{v}(t))-2k_{v}(t-1)\end{aligned},$ (35)

where we have used Eq.(30). Also, based on Eq.(32), it is easy to check the
following equation

$\begin{aligned}
\mathcal{W}_{t}(9)&=\sum_{u,v\in\mathcal{T}_{t-1}}\mathcal{W}_{t}(A_{u}(t);B_{v}(t))+\sum_{v\in\mathcal{T}_{t-1}}\mathcal{W}_{t}(A_{v}(t);B_{v}(t))\\\
&=2\left(\mathcal{W}_{t}(7)-\sum_{v\in\mathcal{T}_{t-1}}W_{t}(B_{v}(t))\right)-2\sum_{u\in\mathcal{T}_{t-1}}\sum_{v(\neq
u)\in\mathcal{T}_{t-1}}k_{v}(t-1)+2\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1)\\\
&=2\mathcal{W}_{t}(7)-2\left(\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1)-2(|\mathcal{T}_{t-1}|-1)\right)-4(|\mathcal{T}_{t-1}|-1)^{2}+2\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1)\\\
&=2\mathcal{W}_{t}(7)-4(|\mathcal{T}_{t-1}|-1)^{2}+4(|\mathcal{T}_{t-1}|-1)\\\
&=8\mathcal{W}_{t}(1)-2(|\mathcal{T}_{t-1}|-1)(4|\mathcal{T}_{t-1}|-1)\end{aligned}.$
(36)

Until now, we have exhaustively enumerated all possible combinatorial ways
when considering distance of an arbitrary pair of vertices in tree
$\mathcal{T}_{t}$. All the unknowns derived in the equations above can be
certainly determined by computing a recurrence that is relevant to vertex
number $|\mathcal{T}_{t-1}|$ and Wiener index $\mathcal{W}_{t-1}$. After that,
putting all things together, we can see that, after performing some simple
arithmetics, Wiener index $\mathcal{W}_{t}$ of scale-free tree network
$\mathcal{T}_{t}$ satisfies

$\mathcal{W}_{t}=\sum_{i=1}^{9}\mathcal{W}_{t}(i)-\mathcal{W}_{t}(3)-\mathcal{W}_{t}(8)=75\mathcal{W}_{t-1}-20|\mathcal{T}_{t-1}|^{2}+20|\mathcal{T}_{t-1}|$
(37)

This means that Wiener index $\mathcal{W}_{t}$ of tree network
$\mathcal{T}_{t}$ is precisely estimated in an iterative manner. And then,
using the elegant relation in Eq.(11), the closed-from solution of mean
hitting time $\langle\mathcal{H}_{t}\rangle$ on tree network $\mathcal{T}_{t}$
is given by

$\langle\mathcal{H}_{t}\rangle=\frac{2}{5^{t}(|\mathcal{T}_{0}|-1)+1}\left\\{75^{t}\mathcal{W}_{0}-4\times
5^{t}\left[\frac{15^{t}-5^{t}}{10}(|\mathcal{T}_{0}|-1)^{2}+\frac{15^{t}-1}{14}(|\mathcal{T}_{0}|-1)\right]\right\\}$
(38)

in which we also make use of Eq.(1).

It is worth mentioning that different from most prior work focusing
particularly on tree networks whose seed is a single edge, formula in Eq.(38)
holds in this sense that an arbitrary tree is selected as the seed for
creating tree network $\mathcal{T}_{t}$. Nonetheless, as an illustrative
example, we plot numerical values for $\langle\mathcal{H}_{t}\rangle$ at
different time steps in Fig.2 where the seed is a single edge. The results
suggest that empirical simulations are in perfect agreement with the
theoretical analysis.

Fig.2. (Color online) The diagram of mean hitting time
$\langle\mathcal{H}_{t}\rangle$ where a single edge serves as the seed of
algorithm $\Phi$. The straight line indicates analytical solutions for
quantity $\langle\mathcal{H}_{t}\rangle$ derived from Eqs.(38). Square
represents the corresponding numerical values produced by computer simulation.
Clearly, theoretical analysis is strongly consistent with experimental
simulation.

In the thermodynamic limit, we can find the next expression

$\langle\mathcal{H}_{t}\rangle\sim|\mathcal{T}_{t}|^{\chi},\quad\chi=1+\frac{\ln
3}{\ln 5}.$ (39)

In the past, such tree models whose mean hitting time exhibits a power-law
function about vertex number are rarely encountered in the study of scale-free
tree networks. Although some fractal trees, for instance, Vicsek fractal [27]
and T-graph [29], show such a phenomenon when studying mean hitting time on
their own underlying structures, it is well known that they have no scale-free
feature. In a nutshell, our tree network $\mathcal{T}_{t}$, as one of
important members in tree family, helps ones to in-depth probe rich underlying
structures of tree models. At the same time, it also provides some instruments
for building some tree networks of great interest.

As mentioned previously, mean hitting time
$\langle\mathcal{H}^{\star}_{t}\rangle$ for random walks on scale-free tree
network $\mathcal{T}^{\star}_{t}$ is exactly derived in a similar manner as
well. In fact, almost all calculations necessary for deriving quantity
$\langle\mathcal{H}^{\star}_{t}\rangle$ have been reported in the process of
computing $\langle\mathcal{H}_{t}\rangle$ and one just needs to make some
slight modifications before reaching his goal. Therefore, we omit the
corresponding detailed derivations and only provide an approximate solution to
quantity $\langle\mathcal{H}^{\star}_{t}\rangle$, i.e.,

$\langle\mathcal{H}^{\star}_{t}\rangle\sim|\mathcal{T}^{\star}_{t}|^{\chi^{\star}}\ln|\mathcal{T}^{\star}_{t}|,\quad\chi^{\star}=1$
(40)

where we have used Eq.(2). From Eqs.(39) and (40), this implies that while
trees $\mathcal{T}_{t}$ and $\mathcal{T}^{\star}_{t}$ have the same vertex
number, their own underlying structures play completely different roles on
studying random walks.

###  4.2 Discussions

Below, we continue to discuss about trees $\mathcal{T}_{t}$ and
$\mathcal{T}^{\star}_{t}$ in more detail for purpose of revealing similarities
and differences. In addition, some other related topics also are provided.
That is to say, this section aims at extending some results in the previous
sections.

First of all, let us rethink of diameter and mean hitting time on two tree
networks $\mathcal{T}_{t}$ and $\mathcal{T}^{\star}_{t}$. In the large graph
size limit, there is an apparent connection, namely,

$\langle\mathcal{H}_{t}\rangle/|\mathcal{T}_{t}|\sim
D_{t},\quad\langle\mathcal{H}^{\star}_{t}\rangle/|\mathcal{T}^{\star}_{t}|\sim
D^{\star}_{t}.$

This tells us that for a tree network $\mathcal{T}$, diameter can be used as
an index for approximately measuring mean hitting time for random walks. One
of important reasons for this is the elegant relationship between Wiener index
$\mathcal{W}_{\mathcal{T}}$ and mean hitting time
$\langle\mathcal{H}_{\mathcal{T}}\rangle$ mentioned in Eq.(11). Roughly
speaking, if one can determine a precise enough solution of diameter of tree
network, then it is not necessary to calculate the corresponding mean hitting
time when he only wants an approximate description about quantity
$\langle\mathcal{H}_{\mathcal{T}}\rangle$.

On the other hand, from the theory point of view, it is interesting yet
challenging to capture exact solution of mean hitting time on various kinds of
trees. It is mainly because in such a rigorous calculational process, one may
develop some effective techniques in order to address these tasks and also
uncover connections some topological parameters behind tree networks. For
instance, the term $\sum_{v\in\mathcal{T}_{t-1}}k^{2}_{v}(t-1)$ in our
calculations is not derived straightforwardly. To this end, we establish
convenient transformation methods from unknowns to some quantities known to
us. As a result, some simple and elegant relationships between distinct
quantities are built, such as, formula in Eq.(37). This further enables
understanding of topological structure of tree networks. In a word, it is
still of great interest to study closed-form solution of some fundamental and
important structural parameters on a verity of tree networks in a
mathematically rigorous manner.

As explained by algorithm $\Phi$ in Section II, tree network $\mathcal{T}_{t}$
is generated via two kinds of generative mechanisms. However, only one of the
both mechanisms is used to create tree network $\mathcal{T}^{\star}_{t}$. Our
analysis shows that doing so leads to two different kinds of trees. Although
they have scale-free feature, the power-law exponents fall into two distinct
regions. As pointed out in our prior work [31], if one only inserts new
vertices on existing edge of tree, the mean hitting time for the end tree can
asymptotically reach to the theoretical upper bound in the large graph size
limit. On the contrary, connecting new vertices to each vertex of tree
directly at each time step leads to a tree whose mean hitting time gradually
tends to the opposite direction, i.e., to the theoretical lower bound. Such an
example is tree $\mathcal{T}^{\star}_{t}$. Now, it is clear to see that the
scaling of mean hitting time on tree network $\mathcal{T}_{t}$ is within the
middle region between theoretical lower bound and theoretical upper bound.
While some fractal trees [27]-[30] possess such a characteristic, they do not
have scale-free feature. On balance, the light shed by constructing our tree
$\mathcal{T}_{t}$ is helpful to generate other interesting trees in the
future.

##  5 conclusion

In summary, we put forward a class of growth tree networks $\mathcal{T}_{t}$
and discuss some fundamental structural parameters. Model $\mathcal{T}_{t}$
turns out to follow power-law degree distribution with exponent larger than
$3$, which is seldom reported in the realm of scale-free tree networks. Next,
small-world property can not be observed on tree $\mathcal{T}_{t}$ because it
has a large diameter, which is completely different from most of published
scale-free trees. Additionally, in order to better understand the underlying
structure of tree $\mathcal{T}_{t}$, we study random walks and then derive an
exact formula for mean hitting time. The results show that as opposed to most
scale-free trees, the mean hitting time of our tree network $\mathcal{T}_{t}$
obeys a power-law function whose variable is vertex number of tree
$\mathcal{T}_{t}$. Note that another scale-free network
$\mathcal{T}^{\star}_{t}$ is in-depth analyzed as well. As mentioned in the
text, tree network $\mathcal{T}_{t}$ shows significantly distinct topological
structure with its analogous $\mathcal{T}^{\star}_{t}$ in terms of several
structural indices while they have a vertex number in common. Last but not the
least, we firmly believe that tree network $\mathcal{T}_{t}$ has many other
intriguing features that are worthy of exploring, which is left as research
direction in our future work.

## Acknowledgments

The research was supported in part by the National Key Research and
Development Plan under grant 2020YFB1805400 and the National Natural Science
Foundation of China under grant No. 61662066.

## References

  * [1]
  * [2] A.-L. Barabási. Network Science. Cambridge University Press. (2016).
  * [3] M.E.J. Newman. The structure and function of complex networks. SIAM Review. 45(3), 167-256 (2003).
  * [4] M.E.J. Newman. Network: An Introduction. Oxford University Press. (2020)
  * [5] D.J. Watts and S.H. Strogatz. Collective dynamics of small-world networks. Nature. 393(6684), 440-442 (1998).
  * [6] A.-L. Barabási and R. Albert. Emergence of scaling in random networks. Science. 286(5439), 509-512 (1999).
  * [7] M. Catanzaro and R. Pastor-Satorras. Analytic solution of a static scale-free network model. The European Physical Journal B. 44, 241-248 (2005).
  * [8] Y.B. Xie, Z. Tao and B.H. Wang. Scale-free networks without growth. Physica A. 387, 1683-1688 (2008).
  * [9] F. Ma, J. Su, Y.X. Hao, B. Yao and G.H. Yan. A class of vertex-edge-growth small-world network models having scale-free, self-similar and hierarchical characters. Physica A. 492, 1194-1205 (2018).
  * [10] L.A.N. Amaral, A. Scala, M. Barthelemy and H.E. Stanley. Classes of Small-World Networks. Proceedings of the National Academy of Sciences. 97, 11149-11152 (2000).
  * [11] Q.K. Telesford, K.E. Joyce, S. Hayasaka, J.H. Burdette and P.J. Laurienti. The ubiquity of small-world networks. Brain Connectivity. 1, 367-375 (2011).
  * [12] R. Cohen and S. Havlin. Scale-Free Networks Are Ultrasmall. Phys. Rev. Lett. 90, 058701 (2003).
  * [13] Z.-Z. Zhang, S.-G. Zhou and T. Zou. Self-similarity, small-world, scale-free scaling, disassortativity, and robustness in hierarchical lattices. The European Physical Journal B. 56, 259-271 (2007).
  * [14] F. Ma, X.M. Wang and P. Wang. An ensemble of random graphs with identical degree distribution. Chaos. 30, 013136 (2020).
  * [15] X.M. Wang and F. Ma. Constructions and properties of a class of random scale-free network. Chaos. 30, 043120 (2020).
  * [16] J.S. Andrade, H.J. Herrmann, R.F.S. Andrade and L.R.D. Silva. Apollonian Networks: Simultaneously Scale-Free, Small World, Euclidean, Space Filling, and with Matching Graphs. Phys. Rev. Lett. 94, 018702 (2005).
  * [17] S.N. Dorogovtsev, A.V. Goltsev and J.F.F Mendes. Pseudofractal scale-free web. Phys. Rev. E. 65, 066122 (2002).
  * [18] E. Ravasz and A.-L. Barabási. Hierarchical organization in complex networks. Phys. Rev. E. 67, 026112 (2003).
  * [19] C.T. Diggans, E.M. Bollt and D. ben-Avraham. Stochastic and mixed flower graphs. Phys. Rev. E. 101, 052315 (2020).
  * [20] S. Athreya, W. Löhr and A. Winter. Invariance principle for variable speed random walks on trees. Annals of Probability. 45, 625-667 (2014).
  * [21] A. Baronchelli, M. Catanzaro and R. Pastor-Satorras. Random walks on complex trees. Phys. Rev. E. 78, 011114 (2008).
  * [22] R.M. D’Souza, P.L. Krapivsky and C. Moore. The power of choice in growing trees. The European Physical Journal B. 59, 535-543 (2007).
  * [23] A. Gabel, P.L. Krapivsky and S. Redner. Highly dispersed networks generated by enhanced redirection. Journal of Statistical Mechanics: Theory and Experiment. 4, P04009 (2014).
  * [24] Z.-M. Lu and S.-Z. Guo. A small-world network derived from the deterministic uniform recursive tree. Physica A. 391, 87-92 (2010).
  * [25] Z.Z. Zhang, Y. Qi, S.G. Zhou, S.Y. Gao and J.H. Guan. Explicit determination of mean first-passage time for random walks on deterministic uniform recursive trees. Phys. Rev. E. 81, 016114 (2010).
  * [26] F. Ma, P. Wang and B. Yao. Random walks on Fibonacci treelike models: emergence of power law. arXiv:1904.11314v2.
  * [27] T. Vicsek. Fractal models for diffusion controlled aggregation. J. Phys. A: Math. Gen. 16, L647-L652 (1983).
  * [28] C.S. Jayanthi, S.Y. Wu and J. Cocks. Real space Green’s function approach to vibrational dynamics of a Vicsek fractal. Phys. Rev. Lett. 69, 1995 (1992).
  * [29] S. Redner. A Guide to First-passage Processes. Cambridge University Press, (2001).
  * [30] E. Agliari. Exact mean first-passage time on the T-graph. Phys. Rev. E. 77, 011128 (2008).
  * [31] F. Ma, P. Wang and X.D. Luo. A Method for Geodesic Distance on Subdivision of Trees With Arbitrary Orders and Their Applications. IEEE Transactions on Knowledge and Data Engineering. DOI: 10.1109/TKDE.2020.3014191 (Early Access Article).
  * [32] M.E. Dyer, A. Galanis, L.A. Goldberg, M.R. Jerrum and E.J. Vigoda. Random Walks on Small World Networks. ACM Transactions on Algorithms. 16, 1-33 (2020).
  * [33] J.D. Noh and H. Rieger. Random walks on complex networks. Phys. Rev. Lett. 92, 118701 (2004).
  * [34] A. Beveridge and M. Wang. Exact Mixing Times for Random Walks on Trees. Graphs and Combinatorics. 29(4), 757-772 (2013).
  * [35] F. Ma and P. Wang. Random walks on a tree with applications. Phys. Rev. E. 102, 022305 (2020).
  * [36] K. Mokhtar, A.A. Fahimah and T. András. Commute times of random walks on trees. Discrete Applied Mathematics. 161(7), 1014-1021 (2013).
  * [37] D.H. Kim, J.D. Noh and H. Jeong. Scale-free trees: the skeletons of complex networks. Phys. Rev. E. 70, 46126 (2004).
  * [38] J.A. Bondy and U.S.R. Murty. Graph Theory. Berlin, Germany: 1170 Springer, 2008.
  * [39] H.K. Dai and M. Toulouse. Lower Bound on Network Diameter for Distributed Function Computation. in Processing of International Conference on Future Data and Security Engineering. 239-251 (2019).
  * [40] G. Bianchin, F. Pasqualetti and S. Zampieri. The role of diameter in the controllability of complex networks. in Processing of 54th IEEE Conference on Decision and Control (CDC). 980-985 (2015).
  * [41] F. Ma, X.M. Wang, P. Wang and X.D. Luo. Dense networks with scale-free feature. Phys. Rev. E. 101, 052317 (2020).
  * [42] C.I. Oliver. Elements of Random Walk and Diffusion Processes. Wiley. (2013).
  * [43] S. Hwang, D.-S. Lee and B. Kahng. First passage time for random walks in heterogeneous networks. Phys. Rev. Lett. 109, 088701 (2012).
  * [44] T.J. Perkins, E. Foxall, L. Glass and R. Edwards. A scaling law for random walks on networks. Nature Communications. 5, 5121 (2014).
  * [45] A.K. Chandra, P. Raghavan, W.L. Ruzzo, R. Smolensky and P. Tiwari. The electrical resistance of a graph captures its commute and cover times. Computational Complexity. 6, 312-340 (1996).
  * [46] A. Georgakopoulos and S. Wagner. Hitting times, cover cost and the Wiener index of a tree. Journal of Graph Theory. 84(3), 311-326 (2017).