Turbulence, driven by density and temperature gradients, plays an important role in strongly magnetized plasmas of magnetic confinement devices, such as tokamaks, as it causes turbulent heat and particle transport.
 However, it also causes the formation of large scale zonal structures, which can in turn reduce the turbulent transport by shearing apart turbulent eddies. In this work, the interplay between turbulence and zonal structures is investigated in gyrokinetic simulations of Ion Temperature Gradient (ITG) driven turbulence, using the Gyrokinetic Electromagnetic Numerical Experiment (GENE).
 By carefully scanning over the ion temperature gradient, the existence of a region of quasi-null heat flux, for temperature gradients exceeding the linear stability threshold, usually referred to as the Dimits non-linear upshift, is recovered.
 In gyrokinetics, the free energy is a quadratic quantity conserved by the non-linear terms: by introducing the zonal free energy fraction, i.e. the ratio between the zonal and total free energy as an order parameter, we show that the Dimits non-linear upshift can be interpreted as a phase transition, between a zonal dominated, ordered, state, and a chaotic, turbulent state.
 Since the control parameter of this transition is the ion temperature gradient, we further show, using a sequence of simulations, that this phase transition can exhibit hysteresis when the temperature gradient is progressively increased or decreased.