How to interpret the output of R's glmer() with quadratic terms

I want to use quadratic terms to fit my general linear mixed model with id as a random effect, using the lme4 package. It’s about how the distance to settlements influences the probability of occurrence of an animal. I use the following code (I hope it is correct):

glmer_dissettl <- glmer(case ~ poly(dist_settlements,2) + (1|id), data=rsf.data, family=binomial(link="logit"))

summary(glmer_dissettl) 

I get the following output:

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: case ~ poly(dist_settlements, 2) + (1 | id)
   Data: rsf.data

     AIC      BIC   logLik deviance df.resid 
  6179.2   6205.0  -3085.6   6171.2     4654 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.14647 -0.90518 -0.04614  0.94833  1.66806 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.02319  0.1523  
Number of obs: 4658, groups:  id, 18

Fixed effects:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 0.02684    0.04905   0.547    0.584    
poly(dist_settlements, 2)1 37.94959    2.41440  15.718   <2e-16 ***
poly(dist_settlements, 2)2 -1.28536    2.28040  -0.564    0.573    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) p(_,2)1
ply(ds_,2)1 0.083         
ply(ds_,2)2 0.067  0.150  

I don’t know exactly how to interpret this, especially with the two lines for poly(dist_settlements,2). Next to understanding, I also wanna see if the quadratic term is making the model better than the basic model without it.

The output of the basic model without a quadratic term:

Generalized linear mixed model fit by maximum likelihood
  (Laplace Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: case ~ scale(dist_settlements) + (1 | id)
   Data: rsf.data

     AIC      BIC   logLik deviance df.resid 
  6177.5   6196.9  -3085.8   6171.5     4655 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6009 -0.8998 -0.0620  0.9539  1.6417 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.02403  0.155   
Number of obs: 4658, groups:  id, 18

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)
(Intercept)              0.02873    0.04945   0.581    0.561
scale(dist_settlements)  0.55936    0.03538  15.810   <2e-16
                           
(Intercept)                
scale(dist_settlements) ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scl(dst_st) 0.077 

I appreciate every tip.

>Solution :

A couple of points.

• Coefficients of non-linear model terms do not have a straightforward interpretation and you should make effect plots to be able to communicate the results from your analyses. You may use effectPlotData() from the GLMMadaptive package to do this. Refer to this page for more information.
• To be able to appraise whether including a quadratic effect of dist_settlements improves the model fit, you should fit a model without the squared term (i.e. only the linear effect of dist_settlements) and a model with the squared term. Make sure to fit both models using maximum likelihood, not REML. Then perform a likelihood ratio test to appraise whether inclusion of complex terms improves the model fit.
• The variance of the random intercepts is rather close to 0, which may require your attention. Refer to this answer and this section of Ben Bolker’s github for more information on this topic.