Thermodynamics of a Higher-Order Topological Insulator

We investigate the order of the topological quantum phase transition in a two dimensional quadrupolar topological insulator within a thermodynamic approach. Using numerical methods, we separate the bulk, edge and corner contributions to the grand potential and detect different phase transitions in the topological phase diagram. The transitions from the quadrupolar to the trivial or to the dipolar phases are well captured by the thermodynamic potential. On the other hand, we have to resort to a grand potential based on the Wannier bands to describe the transition from the trivial to the dipolar phase. The critical exponents and the order of the phase transitions are determined and discussed in the light of the Josephson hyperscaling relation.

We investigate the order of the topological quantum phase transition in a two dimensional quadrupolar topological insulator within a thermodynamic approach. Using numerical methods, we separate the bulk, edge and corner contributions to the grand potential and detect different phase transitions in the topological phase diagram. The transitions from the quadrupolar to the trivial or to the dipolar phases are well captured by the thermodynamic potential. On the other hand, we have to resort to a grand potential based on the Wannier bands to describe the transition from the trivial to the dipolar phase. The critical exponents and the order of the phase transitions are determined and discussed in the light of the Josephson hyperscaling relation.

I. INTRODUCTION
Topological states of matter have attracted a great deal of attention during the last decade. In addition to the usual discrete symmetries that compose the tenfold way [1][2][3], the periodic table of topological materials was extended by the inclusion of crystalline symmetries [4,5], which led to the prediction of new kinds of topological phases. A typical example is the recent proposal of higher-order topological insulators (HOTIs) [6].
HOTIs, or more generally higher-order topological phases, are such that the symmetry protected modes associated with the boundary between the topological phase and the vacuum occur not in the system with codimension 1, but with codimension 2 or higher. This means that for two-dimensional (2D) systems, these modes will appear at the 0D corners. Similarly, for 3D systems they do not arise at the 2D surface but at the 1D hinges or at the 0D corners of the system. The existence of HO-TIs was experimentally verified in a variety of systems, such as bismuth [7], electrical circuits [8], acoustic [9,10] and electronic systems [11], and new types of HOTIs and higher-order topological superconductors were recently proposed [12][13][14][15].
A very interesting question posed by this kind of topological phase is the fate of the bulk-boundary correspondence [16]. A thermodynamic analysis of topological quantum phase transitions (TQPT) [17][18][19][20] inspired by Hill's thermodynamics [21] has been successful in describing the critical behaviour of several models and recently this method was used to devise topological heat machines [22]. Here, we extend this thermodynamical approach to describe HOTIs, i.e., to deal with systems that have both edge and corner modes. For this, we investigate the model proposed in Ref. [6], which is a 2D SSH model, with broken symmetry between the x and y directions and different intra-and inter-cell hopping parameters. In addition, a π-flux pierces the plaquettes, such that some of the hopping parameters are negative. This model exhibits a rich phase diagram, with phases characterised by the polarisation: depending on the ratio between the different hopping parameters, the system is in a quadrupolar, a dipolar or a trivial phase. To describe the phase transitions, we calculate the correlation functions for the bulk, edge and corner contributions, and look for divergences or jumps in their derivatives, which signal a phase transition.
For the transition from the trivial to the quadrupolar state, we find that the bulk exhibits a third-order, the edge a second-order, and the corner a first-order phase transition. The transition from the dipolar to the quadrupolar phase is also visible in the thermodynamic response at the edge and corner contributions, although it is invisible for the bulk. Finally, the transition from the trivial to dipolar phase cannot be captured by the thermodynamic potential description, since there is no band gap closing at the transition. To characterise the topology of these phases, the Wannier spectrum was investigated [6]. The gap in the spectrum of the Wannier centers closes for all phase transitions, and therefore we extracted an effective grand potential from the Wannier spectrum. We found that the Wannier grand potential exhibits discontinuities at all phase transitions, fully characterising the critical behaviour of the system. The critical exponents were extracted by analysing the functional dependence of the closing of either the band or Wannier gap, and the order of the phase transitions were determined on the basis of the Josephson hyperscaling relation. An excellent agreement is found with the order obtained by analysing the discontinuities in the thermodynamic or the Wannier grand potential.
The outline of the paper is as follows: In Section II, we present the quadrupolar topological insulator model proposed in Ref. [6], focusing specifically on the different phases of the phase diagram. Next, in Section III, we extend the thermodynamical description of topolog-ical materials to deal with systems composed of bulk, edge and corner contributions, and explain how one can devise a subtraction scheme to identify the critical behavior of each part. In Section IV, we apply this approach to investigate the quadrupolar to trivial phase transition and show how the critical exponents can be understood by the Josephson's hyperscaling relation. In Section V, we perform the same analysis to other phase transitions in the phase diagram and show that the procedure does not apply due to the absence of a gap closing in some of the topological phase transitions. In Section VI, we show that a grand potential calculated from the Wilson loop spectra, that we name Wannier grand potential, is sensitive to all phase transitions, but with critical exponents of a system with smaller dimension, which we identify to be the edge of the original lattice.

II. QUADRUPOLAR TOPOLOGICAL INSULATOR
The existence of HOTIs was proposed based on the modern theory of polarisation [6,[23][24][25], which associates the projected position operator to the Wilson loop and a topological index. In these systems, there are symmetry protected modes that are localised not on the entire surface of the material, but on parts of this surface. These modes are analogous to multipoles in classical electrodynamics. To realise a quadrupolar topological insulator in a 2D lattice, such that the protected modes are localised in the corners of the lattice, Benalcazar et al. [6] proposed a model which, in presence of inversion symmetry, has a vanishing dipolar order and a nonzero charge localised in the corner due to the bulk quadrupolar moment.
The Hamiltonian of this model reads where c ( †) R,i denotes the destruction (creation) operator at site i in the unit cell labelled by R (see sketch in the inset of Fig. 1), and γ x,y (λ x,y ) are the intracell (intercell) hopping parameters in the x or y directions, respectively. Notice the negative hoppings terms in Eq. (1), which correspond to a π-flux (gauge field) piercing the plaquettes. The model exhibits a quadrupolar phase transition when the hopping parameters are γ x = γ y = γ, λ x = λ y = 1 and |γ| < 1, such that the different corners are related by inversion symmetry and the π-flux on the hopping cancels the edge currents.
The phase diagram of this system is shown in Fig. 1 [6]. The topological index that characterises the topological x/y is the x/y polarisation of the Wannier band ζ x/y (k x /k y ) localised below the gap. This polarisation is quantized in units of 1/2 and is defined mod 1, such that it is either 0 or 1/2, characterizing a Z 2 ×Z 2 invariant. The trivial phase has p = (0, 0), the dipolar phase has p = (1/2, 0) or p = (0, 1/2) and the quadrupolar phase has p = (1/2, 1/2). As the polarisation is calculated from the Wannier bands related to the bulk spectrum, the appearance of zero modes in the corner is related to a gap closing in the bulk, showing a new kind of bulk-boundary correspondence [16].  (1) showing the topological phases. For |γx/λx| ≤ 1 and |γy/λy| ≤ 1 the system is in the quadrupolar phase. The dipolar phases are the ones where only one of the ratios |γx/λx| or |γy/λy| is smaller than 1; when both are larger than 1, the system is in the trivial state. Inset: schematic image of the lattice for Hamiltonian (1). The intercell hopping parameters λ are represented by dashed lines, while the intracell hopping parameters γ are represented by solid lines. Negative values of hopping, created by a π-flux, are represented by red links. We will consider in the following sessions the transitions between: the trivial and quadrupolar state (path I in red), the dipolar to quadrupolar state (path II in blue), and the trivial to dipolar state (path III in green).

III. THERMODYNAMICS OF THE QUADRUPOLAR PHASE TRANSITION
TQPTs are QPTs and therefore one should in principle be able to describe their thermodynamic behaviour. However, the Ginzburg-Landau scheme fails to characterise them because they do not have a local order pa-rameter. In addition, taking the thermodynamic limit is also a problem because one is interested in the edges of the system and these do not contribute in this limit anymore. Nevertheless, the Ehrenfest's classification still holds and allows the characterisation of the phase transition by the order of discontinuity of the grand potential.
Since there is a bulk-boundary correspondence in TQPTs, the thermodynamic description is more subtle because the bulk and the boundary scale differently with the size of the system. This leads to a non-extensive grand potential, as only the part of the grand potential related to the bulk scales with the number of particles in the system.
The solution to this conundrum is provided by an approach inspired by Hill thermodynamics [17][18][19][20][21][22]. Within this method, the system is divided into many subsets and the thermodynamic identity is derived without the assumption of extensivity of the energy. Hence, one can describe finite size system, for which the bulk, the edge and the corner contributions scale differently with the system size. For the 2D system under investigation, the bulk scales with L 2 , the edge with L 1 and the corner with L 0 . The grand potential then can be written as where the ω's are the intensive contributions for each component of the system. The scaling of the different ω shown in the Appendix confirms that indeed they are intensive quantities. x-periodic and (d) y-periodic boundary conditions. We fix λx = λy = 1 and vary the intracell hopping parameters as γx = γy = γ (path I in Fig. 1) for a lattice with 40x40 unit cells. The TQPT occurs for |γ| = 1, at which there is a gap closing and the emergence/disappearance of zero-energy modes. The zero modes are present only with open boundary conditions, confirming that indeed these modes are due to corner states. The similarity between the spectrum of the x-periodic and y-periodic is due to the inversion symmetry of Hamiltonian (1) for the chosen parameter values.
For the Hamiltonian (1), we can write Ω as withN denoting the number operator, µ the chemical potential and β = 1/k B T . One can then perform the calculation for different system sizes and separate the components in Eq. (2), by verifying how Ω scales with the system size [17,18]. However, there is another scheme to obtain these contributions using a more geometrical approach [20]. Depending on the boundary conditions of the system, some of these terms are not present, as suggested from the different spectra in Fig. 2. For instance, if we use periodic boundary conditions, there is no boundary between the system and a trivial vacuum. Hence, there is neither a corner nor and edge contribution to the grand potential. Similarly, if we use periodic boundary conditions only along the x or y direction, the system will have no corner contribution, but just one of the edge contributions. Therefore, we can isolate the different terms in Eq. (2) by subtracting the grand potential of systems with different boundary conditions. We represent pictorially the process of subtraction in Fig. 3.

A. Correlation functions in thermodynamics
Instead of dealing with the grand potential, we will consider its derivative. The derivative is related to correlation functions, which will allow us to identify the phase transition. As a matter of fact, if the control parameter of the phase transition is η, the derivative of Ω with respect to η in Eq. (3) is and this is equal to a correlation function for a Hamiltonian in second quantization. Hence, the discontinuity in Ω is encoded in this correlation function [20]. From the set of eigenstates, we can obtain the correlation functions using a matrix S that relates a vector composed of the eigenstates Ψ and a vector composed of the annihilation operators C and we can write a correlation function for the c operators as where ρ and σ are labels for the c operators and m and n for the energy states. As ψ † is the creation operator of the eigenstates of the Hamiltonian, the correlation function ψ † m ψ n (T, µ) will be simply the occupation of each energy level, given by the Fermi-Dirac distribution multiplied by δ m,n Hence, by knowing the spectrum of the system, one can calculate the correlation functions. In particular, for the quantum phase transition under consideration, we evaluate the correlation function for T = 0 and µ = 0. In this case, the Fermi-Dirac distribution is a step function, which vanishes for ε > ε F = 0 and Eq (7) reduces to In this way, the correlation functions can be obtained by diagonalizingĤ. Notice that for different boundary conditions, the spectrum of the system and the wavefunctions can be different, which ultimately lead to different correlation functions.

B. Correlation functions for the quadrupolar transition
To describe the thermodynamic response of a system divided in bulk, edge and corner, we take the derivatives of each term in Eq. (2). For the Hamiltonian (1), we write the derivative of Ω with respect to γ x = γ y = γ given by Eq. (4) as To obtain the terms of Eq. (2), we identify the different contributions for each boundary condition following the scheme of Fig. 3. Since the system with fully-periodic boundary conditions has no boundary, we assume that the contributions coming from this system contains only Ω bulk , such that where the upper index "per" denotes that the correlation function was calculated using periodic boundary conditions. As the system is periodic, the above correlation function does not depend on the position and the sum will be just L 2 times the correlation function for some point that gives us ω bulk . Then, if we subtract the correlation functions associated with the fully-periodic boundary conditions system from the one with open boundary conditions, we obtain only the edge and corner contributions to the grand potential. Next, we assume that the x-periodic system has only the y-edge and the bulk contributions, and similarly for the y-periodic system, Therefore, the edge term corresponding to the sum of Ω x-edge + Ω y-edge will be Finally, the corner contribution can be obtained when we subtract the edge and bulk parts, Calculating then the correlation functions from Eq. (9) for the open, fully-periodic, x-periodic and y-periodic boundary functions, we can obtain the different contributions in Eq. (2). To sum up, using Eq. (8) with ρ and σ taking the lattice indices of Eq. (9), one can obtain the derivatives with respect to γ of the different contributions to the grand potential in Eq. (2), using the correlations functions with the boundary conditions described in Eqs. (10), (14) and (15), which follow the subtraction scheme of Fig. 3.

IV. TRIVIAL TO QUADRUPOLAR PHASE TRANSITION
We start by analysing the transition from the trivial to the quadrupolar phase, path I in Fig. 1. The derivatives of the different contributions to the grand potential are shown in Fig. 4. For these simulations, we take the intercell hopping parameters to be λ x = λ y = 1 in a lattice of 40 × 40 unit cells and we use interpolation of the correlation functions, such that the higher derivatives can be calculated.
In Fig. 4(a), we see that the first derivative of ω diverges for the corner contribution, while to spot the divergences of the grand potential for the edge and bulk we have to go to the second ( Fig. 4(b)) and third derivative (Fig. 4(c)), respectively, which correspond to the first and second derivative of the correlation function. These results are supported by a finite-size scaling shown in the Appendix. This indicates that we can use the thermodynamic approach to identify this TQPT. In addition, we find that the bulk exhibits a third-order phase transition, the edge a second-order, and the corner a first-order one. We can further understand this behavior by analysing the critical exponents using Josephson's hyperscaling relation.

A. Josephson hyperscaling relation
For quantum phase transitions, the canonical critical exponent α, which is related to the order of the phase transition, is determined by the Josephson hyperscaling relation [26] 2 − α = ν (d + z) , where d is the dimension of the system and the critical exponents z and ν are determined by the way how the band gap ∆ closes, with t = (γ c − γ) /γ c the reduced phase transition parameter and p = k − K the momentum relative to the point K, at which the gap closes.
Then, the Josephson hyperscaling relation implies that As 2 − α denotes how the grand potential scales with γ at the phase transition, this relation predicts that the transition should be third-order for the bulk (d = 2), as we obtained from the correlation function. Due to the bulk-boundary correspondence, the phase transition associated with the gap closing in the bulk is the same that is related to the occurrence of zero modes, which implies that the different contributions to the system have the same z and ν, but different values of d as they have different codimensions: the bulk is 2D, the edge is 1D and the corner 0D. Using these dimensions on Eq. (21), we find that the corner has a first-order phase transition and the edge a second-order phase transition, as we obtained from the correlation-function method.

V. FIXED γx PHASE TRANSITIONS
We now apply the same method to investigate two other phase transitions that can be tuned in this system: the one between the dipolar and the quadrupolar phase (obtained by fixing γ x = 0.5 and varying γ y , path II in Fig. 1) and the transition between the trivial and the dipolar phase (obtained by fixing γ x = 1.5 and varying γ y , path III in Fig. 1).
For fixed γ x , the parameter that tunes the phase transition is γ y . Thus, instead of deriving the Hamiltonian (1) with respect to both γ x and γ y like in Eq. (9), we derive it only with respect to γ y We can obtain the different contributions by replacing this correlation function in Eqs. (10), (14) and (15). In  To better understand why there is no response in the grand potential for these transitions, we analyse how the spectrum of the system changes along these phase transitions. In Fig. 6, we display the spectrum for γ x = 0.5 (a-d) and γ x = 1.5 (e-h), while we vary γ y . For the dipolar to quadrupolar transition, there is a gap closing only at the x-edge (Fig. 6(d)) but not on the y-edge (Fig. 6(c)), whereas the bulk (which has periodic boundary conditions on both directions) does not present a gap closing (Fig 6(b)). This seems to explain why there is no signature of phase transition at the bulk while the edge and corner contributions show discontinuities with nearly half the amplitude of the ones in the trivial to quadrupolar phase transition. For γ x = 1.5, none of the spectra shows a gap closing (Fig. 6(e)-(f)).Although there is no gap closing, the wavefunctions change across the phases (see Fig. 7, where we plot the wavefunctions of the highest energy level in the valence band for each phase). We see that for the dipolar phase with γ x = 1.5 and γ y = 0.5, the wavefunction for the x-periodic boundary conditions is identical to the one for the same boundary condition in the quadrupolar phase, so its localised in one of the y-edges and delocalised along x, see Fig. 7. However, for y-periodic boundary conditions, the wavefunction is similar to the one in the trivial phase, which is delocalised, see Fig. 7. These features are captured in the polarisation.
The polarisation invariant changes across this phase transition without a gap closing due to the breaking of inversion symmetry. This happens because the polarisation is associated to a nested Wilson loop [6]. However, the grand potential is not sensitive to the change in this open x-periodic y-periodic Quad. Phase invariant, as it is not sensitive to the presence or absence of inversion symmetry.
The same occurs for the phase transition from the trivial to the dipolar phase. In this case, the corner, edge and bulk contributions are all flat (green lines in Fig. 5) and do not exhibit any sign of a phase transition. This behavior should be contrasted to the red curves in Fig. 5, which perfectly describe the TQPT from the trivial to the quadrupolar phases.
This observation can be clarified in the context of the Josephson hyperscaling relation: it is impossible to define the critical exponents z and ν because there is no actual gap closing at this phase transition. This has the quite striking implication that we are confronted with a TQPT that does not satisfy even Ehrenfest's classification, with all thermodynamical quantities varying smoothly across the phase transition. This kind of behavior is also seen in other types of phase transitions, like the Berezinskii-Kosterlitz-Thouless [27] one, which has an infinite order.  Fig. 1) for (a) fullyperiodic boundary conditions and (b) x-periodic boundary conditions; for the dipolar to quadrupolar phase transition (path II in Fig. 1) for (c) fully-periodic boundary conditions and (d) x-periodic boundary conditions; and for the trivial to dipolar phase transition (path III in Fig. 1) for (e) fullyperiodic boundary conditions and (f) x-periodic boundary conditions. We see a signature of the phase transition in all the plots. We use a system with 100 unit cells on the nonperiodic direction and 100 k-points to calculate the Wilson loops. Note that in Fig. (a-d) we choose a different gauge in relation to the one used in Fig. (e-d) to make more evident the occurrence of a phase transition.

VI. WANNIER GRAND POTENTIAL
To extract thermodynamical information for this type of phase transitions, one should look at some quantities that do change across the different phase transitions. These are the Wannier centers, eigenvalues of the large Wilson loop. Consider the eigenstates |ψ m (k) with band index m. There is a matrix G (k) associated with these eigenstates, which has the matrix elements where m and n are occupied bands and |∆k| = 2π/L. G is called the Wilson loop element [28]. We can then define the large Wilson loop as the successive multiplication of these elements where C denotes a closed path in the first Brillouin zone that goes from k to k. W is a unitary matrix, such that its eigenvalues are of the form e ı2πζ ; due to topology [6,25], these eigenvalues always come in complex conjugated pairs. As a consequence, the Wannier centers ζ, always come in pairs of positive and negative values. The Wilson loop can be written in terms of an effective Hamiltonian, the so called Wannier Hamiltonian H WC such that its eigenvalues are the Wannier centers (at least mod 1) and yield the Wannier spectrum. There is a gauge degree of freedom in the definition of this Hamiltonian, as we can add a term proportional to the identity to it without changing the Wilson loop. This term can be interpreted as an effective chemical potential. One can, in principle, use this effective Hamiltonian to obtain these Wannier bands and derive an effective grand potential from it, the Wannier grand potential Φ W . To test whether one can see a gap closing or any signature of the phase transitions analysed here, we show in Fig. 8 the Wannier spectrum for this transition both for the fullyperiodic and x-periodic Hamiltonian, considering a path C that goes from k x = −π to k x = π. There are clear signatures in the Wannier spectrum of all the phase transitions, either a gap closing (or appearance of "horns" in the trivial to quadrupolar phase transition) for the fully-periodic boundary condition, or the appearance of in-gap states, which are related to a nontrivial polarization, for x-periodic boundary conditions. Since the exact form of the Wannier Hamiltonian in terms of operators is unknown, we are unable to identify which correlation functions are related to the derivative of Φ W . However, as argued in Refs. [6,29], the Wannier Hamiltonian of Eq. (25) is adiabatically connected to the edge Hamiltonian, such that it should correspond to the thermodynamics of a system with a 1D bulk and a 0D boundary. In this case, we assume a form for Φ W that is composed of a bulk φ bulk and edge φ edge contribution as Wannier Free Energy such that we can obtain both contributions by and For free fermions at zero temperature, the grand potential defined in Eq. (3) is given by where m are all available states.
In an analogous way, the Wannier grand potential is defined by where now the chemical potential µ should be set to 0.5 for the trivial to quadrupolar or dipolar to quadrupolar phase transition, as is evident in the spectra of Fig. 8. We calculate these quantities using the Wannier spectra and the results are presented in Fig. 9. The phase transitions for the 1D bulk are second order, while for the 0D edge these transitions are first order. It is interesting to notice that even in the absence of an actual gap closing in the trivial to quadrupolar phase transition, the Wannier grand potential is sensitive to the phase transition.
As the Wannier gap closes with the same critical exponents as the energy gap of Eq. (17), these two transitions are on the same universality class. We see that in this case, these phase transitions indeed satisfy Josephson's hyperscaling relation given in Eq. (21) for a 1D bulk and a 0D edge. The hyperscaling relation together with the scaling of the Wannier grand potential reinforces the identification of the Wannier Hamiltonian with the edge Hamiltonian.
The fact that we can identify all the phase transitions using the Wannier grand potential indicates that the Wilson loop carries topological information that is not contained in the usual band structure, which is already suggested by the absence of any gap closing or appearance of in-gap states in the spectra of Fig. 6. However, assuming that the Wannier centers are indeed adiabatically connected to the edge spectrum, it is quite surprising that they yield a signature of the phase transition, while the actual spectrum, that contains the edge levels, does not provide this information. A possible answer to this question is that, as pointed out in Ref. [29], the density matrix obtained from the Wannier Hamiltonian is basically a purification of the density matrix obtained by the physical Hamiltonian. In this way, it might be possible to obtain the Wannier grand potential after a partial trace in the grand potential over the degrees of freedom related to the periodic direction that compose the large Wilson loop, but this is beyond the scope of this work.

VII. CONCLUSIONS
We extended the formalism used in Refs. [17][18][19][20]22] to investigate HOTIs, where the bulk-boundary correspondence relates the closure of the band gap with the zero modes that occur at the corners of the system. We numerically calculated the spectrum of Hamiltonian (1), which was proposed in Ref. [6], to identify the discontinuities in the derivative of the grand potential and elucidate the order of the TQPT. We found that for the trivial to quadrupolar phase transition, the bulk exhibits a thirdorder, the edge a second-order and the corner a first-order phase transition, in agreement with Josephson's hyperscaling relation.
However, some of the TQPTs are not captured by the thermodynamic potential description, since those showed no discontinuities at transitions between states with different polarisations. This happens due to the absence of a band gap closing at these phase transitions. We then constructed a Wannier grand potential using the Wannier centers obtained via the Wilson loop, as the Wannier bands are expected to track the TQPT in situations when there is no energy gap closing.
The Wannier grand potential presents critical exponents that have, according to Josephson hyperscaling relation, the same universality class of a system with a 1D bulk and a 0D edge. This, together with the scaling of the Wannier grand potential, confirms the identification of the Wannier Hamiltonian with the edge Hamiltonian.
A deeper understanding of this Wannier grand potential and its implication for evaluations of the entanglement entropy or Otto heat cycles, as well as a more complete comprehension of how one can relate the actual grand potential to the Wannier grand potential are interesting topics for future work.

APPENDIX: FINITE SIZE EFFECTS
Here we comment on the main problem that arises due to finite-size effects, that is the hybridisation between corner modes when the system size is too small and the gap does not close at |γ| = 1. In this case, the corner modes hybridize and move away from zero, while the bulk gap remains open. As the correlation function is calculated from the spectrum, we can use the gap closing as a guideline to verify whether the system is large enough to reveal the phase transition. In Fig. 10, we show the finite-size effects for the first derivative of the bulk, edge and corner contribution for different systems sizes. Although a small kink appears close to γ = 1 for the bulk and edge first derivative, it becomes smaller and broader when we increase the system size, indicating a finite-size effect. In comparison, the divergence of the corner contribution becomes narrower as we increase the system size, indicating that is indeed related to a discontinuity.