  It is known that the loss of phase coherence of Cooper pairs in
two-dimensional (2D) superconductivity corresponds to the unbinding of
vortex-antivortex pairs with the quasi-long-range order (quasi-LRO) in the
order-parameter phase field, described by the Berezinskii-Kosterlizt-Thouless
(BKT) transition of a 2D XY model. Here we show that the second-order Josephson
coupling can induce an exotic superconducting phase in a bilayer system. By
using tensor-network methods, the partition function of the 2D classical model
is expressed as a product of 1D quantum transfer operator, whose eigen-equation
can be solved by an algorithm of matrix product states rigorously. From the
singularity shown by the entanglement entropy of the 1D quantum analogue,
various phase transitions can be accurately determined. Below the BKT phase
transition, an inter-layer Ising long-range order is established at
$T_{Ising}$, and the phase coherence of both intra-layers and inter-layers is
locked together. For two identical layers, the Ising transition coincides with
the BKT transition at a multi-critical point. For two inequivalent layers,
however, there emerges an intermediate quasi-LRO phase ($T_{Ising}<T<T_{BKT}$),
where the vortex-antivortex bindings occur in the layer with the larger
intra-layer coupling, but only half-vortex pairs with topological strings exist
in the other layer, corresponding to the phase coherence of pairs of Cooper
pairs. So our study provides a promising way to realize the charge-4e
superconductivity in a bilayer system.
