Subsystem symmetry-protected topological (SSPT) order is a type of quantum order that is protected by symmetries acting on lower-dimensional subsystems of the entire system. In this paper, we show how SSPT order can be characterized and detected by a constant correction to the entanglement area law, similar to the topological entanglement entropy. Focusing on the paradigmatic two-dimensional cluster phase as an example, we use tensor network methods to give an analytic argument that almost all states in the phase exhibit the same correction to the area law, such that this correction may be used to reliably detect the SSPT order of the cluster phase. Based on this idea, we formulate a numerical method that uses tensor networks to extract this correction from ground state wave functions. We use this method to study the fate of the SSPT order of the cluster state under various external fields and interactions, and find that the correction persists unless a phase transition is crossed, or the subsystem symmetry is explicitly broken. Surprisingly, these results uncover that the SSPT order of the cluster state persists beyond the cluster phase, thanks to a new type of subsystem time-reversal symmetry. Finally, we discuss the correction to the area law found in 3D cluster states on different lattices, indicating rich behaviour for general subsystem symmetries.The modern perspective of quantum phases of matter is based on global patterns of entanglement [1][2][3]. Two quantum states are said to lie in distinct (symmetryprotected) topological phases when they cannot be connected by (symmetric) local unitary evolution. As such, topological phases of matter can be characterized and detected by entanglement-based quantities. A well-known example is the topological entanglement entropy (TEE): For states with non-trivial topological order, the scaling of entanglement entropy exhibits a correction to the area law [4][5][6][7]. That is, for a subregion A of a lattice, the entropy for ground states of gapped, local Hamiltonians takes the general form,for some constants a, γ, where ∂A is the boundary of region A and the dots indicate terms that decay exponentially with |A|. The correction γ takes a uniform non-zero value within non-trivial topological phases as a consequence of the global entanglement patterns. It can therefore be used to detect and characterize topological order and topological phase transitions analytically, numerically, and potentially even experimentally [8][9][10][11][12][13][14][15][16][17].Recently, however, it has been observed that γ may deviate from the expected value due to the presence of symmetry-protected topological (SPT) order localized around ∂A [18-23]. One setting in which this occurs is for states with subsystem SPT (SSPT) order [21,[24][25][26]. Such order is non-trivial only in the presence of subsystem symmetries, which are defined as symmetries that act on rigid lower-dimensional subsystems of the entire system. In cases where these symmetries act on 1D lines spanning a 2D lattice, one may ...